Match.h

Go to the documentation of this file.

#ifndef MATCH_H
#define MATCH_H

#include <string>
#include <vector>
#include "Range.h"
#include <iostream>
#include <list>
#include <set>
#include "bioseq.h"

using namespace std;

// Match.h file

class Match {
   public:
      Match() : q(), t() {}
      Match(const Match &m) : q(m.q), t(m.t) { }
      Match(int qb, int qe, int tb, int te) : q(qb,qe), t(tb,te) {}
      ~Match() {}
      Match& operator=(const Match &m);
      bool operator<(const Match &m) const;
      bool operator==(const Match &m) const { return q==m.q && t==m.t; }


      const Range& queryRange() const { return q; }
      const Range& targetRange() const { return t; }
      int targetBegin() const { return t.begin(); }
      int targetEnd() const { return t.end(); }
      int queryBegin() const { return q.begin(); }
      int queryEnd() const { return q.end(); }
      int queryLength() const { return q.length(); }
      int targetLength() const { return t.length(); }
      char targetDirection() const { return t.direction(); }

      bool onPlus() const { return t.direction() == '+' && q.direction() == '+'; }
      bool onMinus() const { return q.direction() == '-' && t.direction() == '-'; }
      /* the protein is always in the + direction.
       * the DNA can be + or -
       * The following methods will probe the direction of the match on DNA.
       */
      bool targetOnPlus() const { return t.direction() == '+'; }
      bool targetOnMinus() const { return t.direction() == '-'; }
      bool isNull() const { return q.isNull() && t.isNull(); }

      
      int qoverlap(const Match &m, bool sameDirection=true) const {
         return q.overlap(m.q, sameDirection); }
      int toverlap(const Match &m, bool sameDirection=true) const {
         return t.overlap(m.t, sameDirection); }
      pair<int, int> overlap(const Match &m, bool sd=true) const
      { return make_pair(q.overlap(m.q, sd), t.overlap(m.t, sd)); }

      pair<int, int> overlay(const Match &m) const
      { return make_pair(q.overlay(m.q), t.overlay(m.t)); }
      pair<float,float> overlayFraction(const Match &m) const {
         return make_pair(q.overlayFraction(m.q), t.overlayFraction(m.t)); }

      /* the overlap/max(L1,L2) where L1 is the range1's length
       * L2 is range2's length
       */
      //pair<float, float> overlapFraction(const Match &m, bool sameDirection=true) const;
      
      pair<float, float> queryOverlapFraction(const Match &m, bool sameDirection=true) const;
      float queryOverlapFractionMax(const Match &m, bool sameDirection=true) const { pair<float,float> f=queryOverlapFraction(m,sameDirection); 
         return f.first > f.second? f.first : f.second;
      }
      float queryOverlapFractionMerge(const Match &m, bool sameDirection=true) const { return q.overlapFraction(m.q, sameDirection); }
      // return the target overlap fractions relative to this and m object
      pair<float, float> targetOverlapFraction(const Match &m, bool sameDirection=true) const;
      float targetOverlapFractionMax(const Match &m, bool sameDirection=true) const {
         pair<float,float> f=targetOverlapFraction(m,sameDirection);
         return f.first>f.second? f.first : f.second; 
      }
      float targetOverlapFractionMerge(const Match &m, bool sameDirection=true) const {
         return t.overlapFraction(m.t, sameDirection); 
      }

      bool contain(const Match &m, bool sameDirection=true) const {
         return q.contain(m.q, sameDirection) && t.contain(m.t, sameDirection);
      }
      bool qcontain(const Match &m, bool sd=true) const {
         return q.contain(m.q,sd);
      }
      bool tcontain(const Match &m, bool sd=true) const {
         return t.contain(m.t,sd);
      }
      bool similar(const Match &m, bool sd=true) const {
         return q.similar(m.q,sd) && t.similar(m.t,sd); }
      bool coveredBy(const Match &m, bool sd=true, float fraction=0.8) const;

      bool plusContain(const Match &m) const { 
         return q.plusContain(m.q) && t.plusContain(m.t); }
      bool qplustminusContain(const Match &m) const { 
         return q.plusContain(m.q) && t.minusContain(m.t); }

      int queryDistance(const Match &m) { 
         return q.distance(m.q); }
      int targetDistance(const Match &m) { 
         return t.distance(m.t); }


      pair<int,int> combine(const Match &m) {
         return make_pair(q.combine(m.q), t.combine(m.t));
      }
      void setNull() { q=Range(); t=Range(); }
      void setTargetRange(int b, int e) { t.set(b,e); }
      void set(int qb, int qe, int tb, int te) {
         q.set(qb,qe); t.set(tb, te); }


      friend ostream& operator<<(ostream &ous, const Match &m) {
         ous << m.q << " x " << m.t; return ous; }
      ostream& output(ostream &ous, const char delimiter[]="\t") const {
         q.output(ous,delimiter) << delimiter;
         t.output(ous,delimiter);
         return ous;
      }
      string asDelimitedString(const char delimiter[]=",") const {
         return q.asDelimitedString(delimiter) + delimiter
            + t.asDelimitedString(delimiter); }

   protected:
      Range q, t;
};

class M8Match : public Match {
   public:
      M8Match() : Match(), iden(0), len(0), mis(0), gap(0), e(9999), _score(0) {}
      M8Match(const M8Match &m) : Match(m), iden(m.iden), len(m.len), 
                                  mis(m.mis), gap(m.gap), e(m.e), _score(m._score) { }
      M8Match(float identity, int alnlen, int mismatches, int gaps, 
            int qbegin, int qend, int tbegin, int tend, double expect, 
            double score)
         : Match(qbegin,qend,tbegin,tend), iden(identity), len(alnlen), 
           mis(mismatches), gap(gaps),e(expect), _score(score) {}
      M8Match& operator=(const M8Match &m);
      float identity() const { return iden; }
      int alnlen() const { return len; }
      double expect() const { return e; }
      double score() const { return _score; }

      friend ostream& operator<<(ostream &ous, M8Match &m);

      ostream& output(ostream &ous, const char delimiter[]="\t") const;
      string asDelimitedString(const char delimiter[]=",") const;
      
      bool contain(const M8Match& m, bool sameDirection=true, float factor=0.8) const;
      void combine(const M8Match &m);

   protected:
      float iden; 
      int len;   
      int mis, gap;
      double e, _score;
};


class cmpMatchTbeginInt {
   public:
      bool operator()(const Match &m, int x) { return m.targetBegin() < x; }
      bool operator()(int x, const Match &m) { return x < m.targetBegin(); }
};
class cmpMatchTbeginIntPtr {
   public:
      // according to SGI, equal_range use comp(*k, value)
      // in template <class ForwardIterator, class T, class StrictWeakOrdering>
      // pair<ForwardIterator, ForwardIterator>
      // equal_range(ForwardIterator first, ForwardIterator last, 
      // const T& value, StrictWeakOrdering comp);
      bool operator()(const Match *m, int x) { return m->targetBegin() < x; }
      // so the following version was not used by the STL
      // but some place in the document, it seemes to require both
      // versions
      bool operator()(int x, const Match *m) { return x < m->targetBegin(); }
};
class cmpMatchTendInt {
   public:
      bool operator()(const Match &m, int x) { return m.targetEnd() < x; }
      bool operator()(int x, const Match &m) { return x < m.targetEnd(); }
};
class cmpMatchTendIntPtr {
   public:
      bool operator()(const Match* m, int x) { return m->targetEnd() < x; }
      bool operator()(int x, const Match* m) { return x < m->targetEnd(); }
};
/*
class MatchTbeginLessThan {
   public:
      MatchTbeginLessThan(int x) : xpoint(x) { }
      bool operator()(const Match &m) { return m.queryBegin()<xpoint; }
   private:
      int xpoint;
};
*/

class cmpQueryPtr {
   public:
   bool operator()(const Match* m1, const Match* m2) const {
      if (m1->queryRange() == m2->queryRange())
         return m1->targetRange() < m2->targetRange();
      else return m1->queryRange() < m2->queryRange(); }
};

class cmpTargetPtr {
   public:
   bool operator()(const Match* m1, const Match* m2) const {
      if (m1->targetRange() == m2->targetRange())
         return m1->queryRange() < m2->queryRange();
      else return m1->targetRange() < m2->targetRange(); }
};
class cmpTargetDirectionalPtr {
   public:
   bool operator()(const Match* m1, const Match* m2) const {
      int cr= m1->targetRange().compareDirectional(m2->targetRange());
      if (cr == -1) return true;
      else if (cr == 1) return false;
      else return m1->queryRange().directionalLess(m2->queryRange());
   }
};

class beforeMatchQuery {
   public:
      beforeMatchQuery(double fr=0.2) : ovlpfr(fr) { }
      bool operator()(const Match &m1, const Match &m2) const {
         return m1.queryRange() < m2.queryRange() && m1.queryOverlapFractionMerge(m2) < ovlpfr;
      }
   private:
      double ovlpfr;
};
class afterMatchQuery {
   public:
      afterMatchQuery(double fr=0.2) : ovlpfr(fr) { }
      bool operator()(const Match &m1, const Match &m2) const {
         return m1.queryRange() > m2.queryRange() && m1.targetOverlapFractionMerge(m2) < ovlpfr;
      }
   private:
      double ovlpfr;
};
class notlessMatchQueryBegin {
   public:
      bool operator()(const Match &m1, const Match &m2) const {
         return m1.queryBegin() >= m2.queryBegin();
      }
};

class M8MatchEX : public M8Match {
   public:
      friend class Boxchain;
      M8MatchEX() : M8Match(), maxScore(0), before(0), after(0) { }
      M8MatchEX(const M8MatchEX &m) : M8Match(m), maxScore(m.maxScore), before(m.before), after(m.after) { }

      /* assign the maxScore to the score of regular blast */
      M8MatchEX(const M8Match &m) : M8Match(m), maxScore(m.score()), before(0), after(0) { }
      M8MatchEX& operator=(const M8MatchEX &m);
      /* maxScore initialized to m._score */
      M8MatchEX& operator=(const M8Match &m);
      double getMaxScore() const { return maxScore; }
      void setMaxScore(double sc) { maxScore=sc; }
      /* only link to previous nodes */
      void attachTo(M8MatchEX &m) {
         maxScore=m.maxScore + score(); before=&m; }
      friend ostream& operator<<(ostream &ous, const M8MatchEX &m);
      M8MatchEX* previous() const { return before; }
      M8MatchEX* next() const { return after; }
      void setPrevious(M8MatchEX* p) { before=p; }
      void setNext(M8MatchEX* p) { after=p; }
      void setDefault() { before=0; maxScore=score(); }


   private:
      double maxScore;
      M8MatchEX* before;
      M8MatchEX* after; // for the chaining algorithm
};
// must be defined after M8MatchEX
class MatchQueryLessThan {
   public:
      MatchQueryLessThan(const M8MatchEX& m, double fr=0.2) 
         : mm(m), ovlpfr(fr) { }
      MatchQueryLessThan(const M8MatchEX* m, double fr=0.2) 
         : mm(*m), ovlpfr(fr) { }
      bool operator()(const M8MatchEX *m) const {
         //return m->queryRange() < mm->queryRange() && m->queryOverlapFractionMerge(*mm) < ovlpfr;
         return m->queryRange() < mm.queryRange() && m->queryOverlapFractionMerge(mm) < ovlpfr;
      }
   private:
      const M8MatchEX& mm;
      double ovlpfr;
};

class ltMatchQueryBeginPtr {
   public:
      bool operator()(const M8MatchEX *p1, const M8MatchEX *p2) const {
         // this version will keep multiple scores for the 
         // same queryBegin location, we don't want this
         /*
         if (p1->queryBegin() == p2->queryBegin()) 
            return p1->getMaxScore() < p1->getMaxScore();
         else return p1->queryBegin() < p2->queryBegin();
         */
         return p1->queryBegin() < p2->queryBegin();
      }
};

class Boxchain;
class M8MatchChain {
   public:
      friend class Tblastn;
      M8MatchChain() : chain(), _first(true), _last(true), sumiden(0), 
                       sumaln(0), sumscore(0) {}
      //M8MatchChain(Boxchain &ch);
      M8MatchChain(M8Match *m) 
         : chain(), _first(true), _last(true), 
           sumiden(m->identity()*m->alnlen()),
           sumaln(m->alnlen()), sumscore(m->score()) { 
            chain.push_back(m); 
            qchain.add(m->queryRange(), true); }

      void add(M8Match *m);
      void nofirst() { _first=false; }
      void nolast() { _last=false; }
      bool hasFirst() const { return _first; }
      bool hasLast() const { return _last; }

      // begining of the model, is on the plus strand, the smallest 
      // if - strand the largest number
      int getTargetBegin() const;
      int getTargetEnd() const;
      int getQueryBegin() const;
      int getQueryEnd() const;

      // this method only make sense if the chain is 
      // not empty, if empty return ' '
      char direction() const;
      int queryOverlap() const { return qchain.getOverlap(); }
      int queryLength() const { return qchain.length().first; }

      // use the M8Match operator <<
      friend ostream& operator<<(ostream &ous, const M8MatchChain &c);
      // not useful for output
      ostream& output(ostream &ous, char delimiter[]="\t") const;

      float avgIdentity() const { return (float)(sumiden/sumaln); }
      double sumScore() const { return sumscore; }
      int size() const { return chain.size(); }

   private:
      vector<M8Match*> chain; // sorted on genomic from small to large

     /* record completness regarding the matches
      * available. Not the true start or end
      */
      bool _first; // contain or overlap first match
      bool _last;  // contain or overlap last match
      RangeChain qchain; // query range chain
      // not needed void computestat();

      double sumiden;
      int sumaln;
      double sumscore;
      //float avgiden;
};

class Tblastn {
   public:
      Tblastn() : qid(), tid(), matches(), matchesQS(), sortedByQuery(false), sortedByTarget(false) {}
      Tblastn(const string &query, const string &target, int queryLen, int targetLen)
         : qid(query), tid(target), qlen(queryLen), tlen(targetLen), 
            matches(), matchesQS(), sortedByQuery(false),
            sortedByTarget(false) {}
      Tblastn(const Tblastn& tbn);
      Tblastn& operator=(const Tblastn& tbn);
      ~Tblastn();
      void clear() { removeAll(); }


      void addMatch(float identity, int alnlen, int mismatches, int gaps,
            int qbegin, int qend, int tbegin, int tend, double expect, double score);
      void addMatch(M8Match *mptr) { matches.push_back(new M8Match(*mptr)); }
      //void removeMatch(M8Match *mptr) { delete mptr; matches.remove(mptr); }
      void add(const Tblastn& tbn);

      void setSortedByTarget(bool yes=true) { sortedByTarget=yes; }
      bool isSortedByTarget() const { return sortedByTarget; }


      // dump all rows to ostream ous
      ostream& output(ostream& ous, char delimiter[]="\t") const;

      void outputModel(const M8MatchChain &amodel, int &id, 
            ostream& ouss, ostream &ousd, char delimiter[]="\t") const;
      void outputModel(const Boxchain &amodel, int &id, 
            ostream& ouss, ostream &ousd, char delimiter[]="\t") const;
      // for human to read
      friend ostream& operator<<(ostream &ous, const Tblastn &b);
      // debug show
      void show(ostream &ous) const;

      
      //int removeJunk(ostream &ous, list<string> &rows, 
       //     const Protein &qseq, const DNA &tseq, float factor=0.8);
      int removeJunk(const Protein &qseq, const DNA &tseq, float factor=0.8);

      void findFootprint(ostream &SUMM, ostream &DETA, float covcut=0.4);
      /* use the sparse dynamic programming to select the best
       * models.  Could be more than one models
       */
      //void bestModel();
      //
      vector<Match> findSections();


      bool empty() { return matches.empty(); }
      const string& query() const { return qid; }
      const string& target() const { return tid; }

      int numMatches() const { return matches.size(); }
      const vector<M8Match*>& getMatches() const { return matches; }
      int queryLength() const { return qlen; }
      int targetLength() const { return tlen; }
      const Range& firstQueryRange() const;
      const Range& lastQueryRange() const;
      float qcov() const;

      bool isFirstQueryMatch(const M8Match* m) const;

      bool overlapFirstQueryMatch(const M8Match* m, int cutoff=9) const;
      bool isLastQueryMatch(const M8Match* m) const;
      bool overlapLastQueryMatch(const M8Match* m, int cutoff=9) const;

      //
      //vector<M8Match*> checkContain() const;
      int checkContain(ostream &ous, float factor=0.8) const;
      bool sameQTPair(const Tblastn& tbn) {
         return qid==tbn.qid && tid == tbn.tid; }


      void sortMatchByTarget();

      void directionalSortMatchByTarget();

      static void readConfig(const string &conf);
      static double getIntronProbility(int len);

   private: 
      string qid, tid; // query and target id
      int qlen, tlen; // query and target length
      vector<M8Match*> matches;
      bool sortedByTarget; // state variable for matches
      void erase(int b, int e);
      int eraseMatch(const Range& r);
      void removeNullMatch();

      mutable vector<M8Match*> matchesQS; // sorted by query
      mutable bool sortedByQuery;
      void sortMatchByQuery() const;
      int removeAll(list<string> &rmlist, ostream &ous);
      int removeAll();
      //vector<M8MatchChain*> models;

      static float qcovCutoff;
      //static vector<pair<int, double> > intronProb;
      static vector<double> intronProb;
};

class Boxchain {
   public:
      //Boxchain() : chain() { }
      Boxchain() : begin(0), end(0), chlen(0), iden(0), qolp(-1) { }
      ~Boxchain();
      Boxchain(M8MatchEX* p);
      Boxchain(const Boxchain &bc);
      //Boxchain& operator=(M8MatchEX* end);
      Boxchain& operator=(const Boxchain& bc);
      void show(ostream &ous) const;
      const M8MatchEX* last() const { return end; }
      const M8MatchEX* first() const { return begin; }
      M8MatchEX* first() { return begin; }
      M8MatchEX* last() { return end; }
      int numMatch() const { return chlen; }
      char direction() const { return begin->targetDirection(); }
      // oputput functions 
      int queryBegin() const { return begin->queryBegin(); }
      int queryEnd() const { return end->queryEnd(); }
      int targetBegin() const { return begin->targetBegin(); }
      int targetEnd() const { return end->targetEnd(); }
      double sumScore() const { return end->getMaxScore(); }
      double maxScore() const { return end->getMaxScore(); }
      double avgIdentity() const;
      int queryOverlap() const;
      ostream& summaryOutput(ostream &ous, const char dl[]="\t") const;
      ostream& detailOutput(ostream &ous, const string &qid, const string &tid,
            const int &id, const char dl[]="\t") const;

      int queryCovered() const { return end->queryEnd()-begin->queryBegin()+1; }
      int querySeqCovered() const;


   private:
      //list<M8MatchEX*> chain;
      M8MatchEX* begin;
      M8MatchEX* end;
      int chlen;  // number of segments
      mutable double iden;
      mutable int qolp;
};

class Linkmatch {
   public:
      ~Linkmatch();
      //TblastnEx(const vector<M8Match*> &match);
      Linkmatch(Tblastn &tb);
      typedef multiset<M8MatchEX*, ltMatchQueryBeginPtr>::iterator bestiterator;
      typedef multiset<M8MatchEX*, ltMatchQueryBeginPtr>::reverse_iterator bestreverse_iterator;

      /* the following two functions reduce performance
       * but they look nicer
       */
      //vector<const M8MatchEX*> startAt(int x) const;
      //vector<M8MatchEX*> endAt(int x) const;
      template<class ForwardIterator> 
         M8MatchEX* bestModel(ForwardIterator first, ForwardIterator last, const Match &box);
      M8MatchEX* bestModel(const Match &sec);
      //void bestModel(vector<Boxchain>& topmodels);
      void bestModel(vector<Boxchain*>& topmodels);  // pointer version
      //void nextBestModel();
      //void bestModel(const Range& sec);
      //
      void showBest(ostream &ous) const;
      void show_best(ostream &ous) const;
      void show_end2match(ostream &ous) const;
      void show_start2match(ostream &ous) const;
      void show_sections(ostream &ous) const;
      //vector<Range> findSections() const;
      void removeChain(M8MatchEX* p, Match &box);

   private:
      //vector<M8MatchEX> elem;
      //list<int> xpoints;
      set<int> xpoints;
      // vector is hard to implement search of Ranges
      multimap<int, M8MatchEX*> start2match;
      multimap<int, M8MatchEX*> end2match;
      multiset<M8MatchEX*, ltMatchQueryBeginPtr> best;
      void cleanUpBest(bestiterator bi);
      vector<Match> sections;
};

#endif