Coverdepth.h

Go to the documentation of this file.

#ifndef COVERDEPTH_H
#define COVERDEPTH_H

#include "GenModel.h"
#include <iostream>
#include <cassert>
#include <limits>

class NegativeCoverage { };

void computeDerivative(const int *in, int *out, int L, int dx=1) 
   throw(NegativeCoverage);
void sumNumber(const int *A, int *&B, int L, int w);

class Dip {
   public:
      Dip() : posTL(-1), posBL(-1), posBR(-1), posTR(-1),
         LH(0), RH(0), BLH(numeric_limits<int>::max()), 
         BRH(numeric_limits<int>::max()) { }
      Dip(int tl, int bl, int br, int tr, int lh, int blh, int brh, int rh) 
         : posTL(tl), posBL(bl), posBR(br), posTR(tr),
           LH(lh), BLH(blh), BRH(brh), RH(rh) { }
      //void setRightDropoff(int pTL, int pL, int pR, int tlH, lH, rH);
      //void setLeftDropoff(int pL, int pR, int pTR, int lH, int rH, int trH);

      friend ostream& operator<<(ostream& ous, const Dip &v);
      void incrementPosition(int B);
      bool isGood() const;
      bool betterThan(const Dip &dip) const;
      int bottomLeft() const { return posBL; }
      int bottomRight() const { return posBR; }
      int bottomWidth() const { return posBR - posBL + 1; }
      int height() const { return max(LH,RH)-min(BLH,BRH); }
      int leftHeight() const { return LH-BLH; }
      int rightHeight() const { return RH-BRH; }
      int bottomHeight() const { return min(BLH,BRH); }
      pair<int,int> breakPoint(char dir='+') const;
      int ratio() const { return max(LH, RH)/min(BLH,BRH); }
      float fineratio() const { return (float)max(LH, RH)/min(BLH,BRH); }
      float fineRatio() const { return (float)max(LH, RH)/min(BLH,BRH); }
      bool isNull() const { return posTL==-1 || posBL==-1 || posBR==-1 || posTR==-1; }
      int posTL, posBL, posBR, posTR;
      int LH, RH; 
      //int LH, RH, BH; 
      int BLH, BRH;
      static int intronAvoid; // 20
};

class RightDropoff : public Dip {
   public:
      RightDropoff() : Dip() {}
      RightDropoff(int pTL, int pL, int pR, int tlH, int lH, int rH) 
         : Dip(pTL, pL, pR, pR, tlH, lH, rH, rH) { }
      int height() const { return LH-BRH; }
      float fineRatio() const { return float(LH)/BRH; }
      int ratio() const { return LH/BRH; }
      bool better(const RightDropoff &drf) { return fineRatio() > drf.fineRatio() && BLH <= drf.BLH; }
      bool isGood() const { return (ratio()>30 && BRH<5) || ratio()>60; }
};

class LeftDropoff : public Dip {
   public:
      LeftDropoff() : Dip() {}
      LeftDropoff(int pL, int pR, int pTR, int blH, int brH, int trH) 
         : Dip(pL, pL, pR, pTR, blH, blH, brH, trH) { }
      int height() const { return RH-BLH; }
      float fineRatio() const { float(RH)/BLH; }
      int ratio() const { return RH/BLH; }
      //int ratio() const;
      bool better(const LeftDropoff &drf) { return fineRatio() > drf.fineRatio() && BRH <= drf.BRH; }
      bool isGood() const { return (ratio()>30 && BLH<6) || ratio()>60; }
      //bool isGood() const;
};

class Coverdepth {
   public:
      Coverdepth() : B(0), E(0), dep(0), direction('?'), der1(0)  { }
      Coverdepth(const Noschain &chain, int C=1);
      Coverdepth(const string &cdstr);
      Coverdepth(int *prof, int bb, int ee, char dir)
         : B(bb), E(ee), dep(prof), der1(0), direction(dir) { }
      ~Coverdepth();
      //~Coverdepth() { delete[] dep; }
      Coverdepth& assign(const Noschain &chain, int C);
      //Coverdepth& operator+(const Noschain &chain);
      Coverdepth& append(const Noschain &chain, int C=1);
      int length() const { if (dep == 0) return 0; return E-B+1; }
      int begin() const { return B; }
      int end() const { return E; }
      bool empty() const { return dep == 0; }
      void clear();
      friend ostream& operator<<(ostream &ous, const Coverdepth &cd);
      Noschain exons() const;
      char strand() const { return direction; }
      void setStrand(char dir) { direction=dir; }
      void reverseStrand() { if (direction=='+') direction='-';
         else if (direction=='-') direction='+'; }
      void trim(const Range &bound);
      void erase(const Range &bound);
      void erase(int bb, int ee);
      void fill(int bb, int ee, int height);
      void fillRight(int bb, int ee);
      void extendRight(int bb, int len);
      void fillLeft(int bb, int ee);
      void extendLeft(int ee, int len);
      bool findDip(int b, int e, Dip &dip) const throw(NegativeCoverage);
      bool searchDip(int b, int e, Dip &dip) const throw(PointOutChain);
      pair<float,float> relativeHeight(const Coverdepth &bcp) const;

   private:
      int findDropoff(list<pair<double, int> > &derivseg, int s, int f, Dip &dip) const;
      bool findRightDropoff(list<pair<double, int> > &derivseg, int s, int f, Dip &dip) const;
      bool findLeftDropoff(list<pair<double, int> > &derivseg, int s, int f, Dip &dip) const;
      //bool findRightDropoff(const vector<pair<double, int> > &derivseg, int s, int f, int &cut) const;
      //int findDropoff(list<pair<double, int> > &derivseg, int s, int f, LeftDropoff &ldrop, RightDropoff &rdrop) const;
      int digDip(int b, int e, vector<Dip> &dips, const int *deriv) const;
      int findVDip(const list<pair<double,int> > &seg, int b, vector<Dip> &dips) const;
      int localMax(int &x, int w=10) const;
      int localMin(int &x, int w=10) const;
      int indexGenomic2dep(int gpos) const { return gpos-B; }
      int indexdep2Genomic(int idx) const { return B+idx; }

   public:
      int* derivative(int dx=1);
      ostream& toPlot(ostream &ous) const;
      float averageHeight() const;
      int area() const;
      int area(int b, int e) const;
      int maxHeight() const;
      pair<int,int> maxAndSum() const;
      void maxAndSum(int &maxh, int &sum) const;
      int summary(int &maxc, int &maxi, int &minc, int &mini, int b=0, int e=-1) const;
      Coverdepth* subprofile(int b, int e, char dir='?') const;
      void truncateTail(int e) { E=e; }
      void truncateHead(int b);

   private:
      bool makeDip(Dip &dip, const int* deriv, const int s, const int ss, const int ff, int &i) const;
      void produceDip(Dip &dip, int dipx[], int b) const;
      void produceVDip(Dip &dip, int dipx[], int b) const;
      int *dep;
      int B, E; // inclusive [B,E]
      char direction; // + or - or ?
      int *der1;
      static int dipedgew;
};

#endif