src/fragments.h

#ifndef FRAGMENTS_H
#define FRAGMENTS_H

#include <string>
#include <vector>
#include <iostream>
#include <fstream>
#include <cassert>
#include <api/BamAlignment.h>
#include "sequence.h"
#include "boost/scoped_ptr.hpp"

typedef size_t TargID;
struct Library;
class Target;
class TargetTable;

enum PairStatus { PAIRED, LEFT_ONLY, RIGHT_ONLY };

struct Indel {
  size_t pos;
  size_t len;
  Indel(size_t p, size_t l) : pos(p), len(l) {}
};

struct ReadHit {
  std::string name;
  bool first;
  bool reversed;
  size_t targ_id;
  size_t left;
  size_t right;
  SequenceFwd seq;
  std::vector<Indel> inserts;
  std::vector<Indel> deletes;
  BamTools::BamAlignment bam;
  std::string sam;
  int mate_l;
};

struct HitParams {
  double align_likelihood;
  double full_likelihood;
  double posterior;
};

class FragHit {
  Target* _target;
  boost::scoped_ptr<ReadHit> _read_l;
  boost::scoped_ptr<ReadHit> _read_r;
  std::vector<const Target*> _neighbors;
  HitParams _params;
  
public:
  FragHit(ReadHit* h) : _target(NULL) {
    if (h->reversed) {
      _read_r.reset(h);
    } else {
      _read_l.reset(h);
    }
  }
  FragHit(ReadHit* l, ReadHit* r) : _target(NULL), _read_l(l), _read_r(r) {
    assert(!l->reversed);
    assert(r->reversed);
    assert(l->name == r->name);
    assert(l->targ_id == r->targ_id);
    assert(l->left <= r->left);
    assert(l->first != r->first);
  }
  std::string frag_name() const {
    if (_read_l) {
      return _read_l->name;
    }
    assert(_read_r);
    return _read_r->name;
  }
  HitParams* params() { return &_params; }
  const HitParams* params() const { return &_params; }
  Target* target() const {
    return _target;
  }
  void target(Target* target) {
    _target = target;
  }
  const std::vector<const Target*>* neighbors() const { return &_neighbors; }
  void neighbors(const std::vector<const Target*>& neighbors) {
    _neighbors = neighbors;
  }
  TargID target_id() const {
    if (_read_l) {
      return _read_l->targ_id;
    }
    assert(_read_r);
    return _read_r->targ_id;
  }
  size_t left() const {
    if (_read_l) {
      return _read_l->left;
    }
    assert(_read_r);
    return _read_r->left;
  }
  size_t right() const {
    if (_read_r) {
      return _read_r->right;
    }
    assert(_read_l);
    return _read_l->right;
  }
  size_t length() const {
    if (_read_l && _read_r) {
      return _read_r->right - _read_l->left;
    }
    return 0;
  }
  const ReadHit* left_read() const {
    if (_read_l) {
      return _read_l.get();
    }
    return NULL;
  }
  const ReadHit* right_read() const {
    if (_read_r) {
      return _read_r.get();
    }
    return NULL;
  }
  const ReadHit* first_read() const {
    if (_read_l && _read_l->first) {
      return _read_l.get();
    }
    assert(_read_r);
    return _read_r.get();
  }
  const ReadHit* second_read() const {
    if (_read_l && !_read_l->first) {
      return _read_l.get();
    } else if (_read_r && !_read_r->first) {
      return _read_r.get();
    } else {
      return NULL;
    }
  }
  ReadHit* left_read() {
    return const_cast<ReadHit*>(const_cast<const FragHit*>(this)->left_read());
  }
  ReadHit* right_read() {
    return const_cast<ReadHit*>(const_cast<const FragHit*>(this)->right_read());
  }
  ReadHit* first_read() {
    return const_cast<ReadHit*>(const_cast<const FragHit*>(this)->first_read());
  }
  ReadHit* second_read() {
    return const_cast<ReadHit*>(const_cast<const FragHit*>(this)->second_read());
  }
  PairStatus pair_status() const {
    if (_read_l && _read_r) {
      return PAIRED;
    }
    if (_read_l) {
      return LEFT_ONLY;
    }
    return RIGHT_ONLY;
  }
};

class Fragment {
  std::vector<FragHit*> _frag_hits;
  std::vector<ReadHit*> _open_mates;
  std::string _name;
  double _mass;
  Library* _lib;
  void add_open_mate(ReadHit* om);

public:
  Fragment(Library* lib);
  ~Fragment();
  const Library* lib() { return _lib; }
  bool add_map_end(ReadHit* r);
  const std::string& name() const { return _name; }
  size_t num_hits() const { return _frag_hits.size(); }
  FragHit* operator[](size_t i) const { return _frag_hits[i]; }
  const std::vector<FragHit*>& hits() const { return _frag_hits; }
  const FragHit* sample_hit() const;
  void mass(double m) { _mass = m; }
  double mass() const { return _mass; }
  void sort_hits();
  bool paired() const {
    if (_frag_hits.empty()) {
      return false;
    }
    return (_frag_hits[0]->pair_status() == PAIRED);
  }
};

#endif