source: branches/0.12-stable/yat/utility/Aligner.cc @ 3260

// $Id: Aligner.cc 3260 2014-06-20 21:09:29Z jari $

/*
  Copyright (C) 2012, 2014 Peter Johansson

  This file is part of the yat library, http://dev.thep.lu.se/yat

  The yat library is free software; you can redistribute it and/or
  modify it under the terms of the GNU General Public License as
  published by the Free Software Foundation; either version 3 of the
  License, or (at your option) any later version.

  The yat library is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
  General Public License for more details.

  You should have received a copy of the GNU General Public License
  along with yat. If not, see <http://www.gnu.org/licenses/>.
*/

#include <config.h>

#include "Aligner.h"
#include "Cigar.h"
#include "Matrix.h"

#include <cassert>
#include <limits>
#include <ostream>

namespace theplu {
namespace yat {
namespace utility {

  Aligner::Aligner(double gap, double open_gap)
    : vertical_gap_(gap), open_vertical_gap_(open_gap), horizon_gap_(gap),
      open_horizon_gap_(open_gap)
  {
  }


  Aligner::Aligner(double vertical_gap, double open_vertical_gap,
                   double horizon_gap, double open_horizon_gap)
    : vertical_gap_(vertical_gap), open_vertical_gap_(open_vertical_gap),
      horizon_gap_(horizon_gap), open_horizon_gap_(open_horizon_gap)
  {
  }


  double Aligner::needleman_wunsch(const Matrix& s)
  {
    Matrix m(s.rows()+1,s.columns()+1, -std::numeric_limits<double>::max());
    m(0, 0) = 0.0;
    // Init upper and left border of matrix
    for (size_t i=1; i<m.rows(); i++)
      m(i,0) = -i * vertical_gap_ - open_vertical_gap_;
    for (size_t i=1; i<m.columns(); i++)
      m(0,i) = -i * horizon_gap_ - open_horizon_gap_;
    (*this)(s, m);
    // return lower right corner element
    return m(s.rows(), s.columns());
  }


  double Aligner::smith_waterman(const Matrix& s)
  {
    Matrix m(s.rows()+1,s.columns()+1);
    (*this)(s, m);
    return max(m);
  }


  void Aligner::operator()(const Matrix& score, Matrix& x)
  {
    assert(x.columns() == score.columns()+1);
    assert(x.rows() == score.rows()+1);

    columns_ = x.columns();
    alignment_.resize(columns_*x.rows(), none);

    for (size_t i=1; i<x.rows(); ++i)
      for (size_t j=1; j<x.columns(); ++j) {
        // match (or mismatch)
        if (x(i-1, j-1) + score(i-1, j-1) > x(i,j)) {
          x(i,j) = x(i-1, j-1) + score(i-1, j-1);
          directions(i, j) = diagonal;
        }

        if (directions(i, j-1) == right) {
          if (x(i, j-1)-horizon_gap_ > x(i,j)) {
            x(i,j) = x(i,j-1) - horizon_gap_;
            directions(i,j) = right;
          }
        }
        else if(x(i,j-1)-horizon_gap_-open_horizon_gap_ > x(i,j)) {
          x(i,j) = x(i,j-1) - horizon_gap_ - open_horizon_gap_;
          directions(i,j) = right;
        }

        if (directions(i-1,j) == down) {
          if (x(i-1,j) - vertical_gap_ > x(i,j)) {
            x(i,j) = x(i-1,j) - vertical_gap_;
            directions(i,j) = down;
          }
        }
        else if (x(i-1,j) - vertical_gap_ - open_vertical_gap_ > x(i,j)) {
          x(i,j) = x(i-1,j) - vertical_gap_ - open_vertical_gap_;
          directions(i,j) = down;
        }
      }
  }


  const Aligner::direction& Aligner::alignment(size_t i, size_t j) const
  {
    return alignment_[i*columns_ + j];
  }


  const Aligner::Cigar Aligner::cigar(size_t row, size_t column) const
  {
    Cigar cig;
    while (row || column) {
      switch (alignment(row, column)) {
      case(diagonal):
        assert(row);
        --row;
        assert(column);
        --column;
        cig.push_front(BAM_CMATCH);
        break;
      case (right):
        assert(column);
        --column;
        cig.push_front(BAM_CINS);
        break;
      case (down):
        assert(row);
        --row;
        cig.push_front(BAM_CDEL);
        break;
      default:
        return cig;
      }
    }
    return cig;
  }


  Aligner::direction& Aligner::directions(size_t i, size_t j)
  {
    return alignment_[i*columns_ + j];
  }


  // ================ Cigar =======================

  void Aligner::Cigar::clear(void)
  {
    cigar_.clear();
  }


  uint8_t Aligner::Cigar::op(size_t i) const
  {
    assert(i<cigar_.size());
    return bam_cigar_op(cigar_[i]);
  }


  uint32_t Aligner::Cigar::length(uint8_t mask) const
  {
    uint32_t res = 0;
    for (size_t i=0; i<cigar_.size(); ++i)
      if (bam_cigar_type(op(i)) & mask)
        res += oplen(i);
    return res;
  }


  uint32_t Aligner::Cigar::length(void) const
  {
    return length(3);
  }


  char Aligner::Cigar::opchr(size_t i) const
  {
    assert(i<cigar_.size());
    return bam_cigar_opchr(cigar_[i]);
  }


  uint32_t Aligner::Cigar::oplen(size_t i) const
  {
    assert(i<cigar_.size());
    return bam_cigar_oplen(cigar_[i]);
  }


  void Aligner::Cigar::pop_back(uint32_t len)
  {
    while (len) {
      assert(size());
      size_t i = size()-1;
      if (len < oplen(i)) {
        cigar_[i] -= (len<<BAM_CIGAR_SHIFT);
        return;
      }
      else {
        len -= oplen(i);
        cigar_.pop_back();
      }
    }
  }


  void Aligner::Cigar::pop_front(uint32_t len)
  {
    while (len) {
      assert(size());
      if (len < oplen(0)) {
        cigar_[0] -= (len<<BAM_CIGAR_SHIFT);
        return;
      }
      else {
        len -= oplen(0);
        cigar_.pop_front();
      }
    }
  }


  void Aligner::Cigar::push_back(uint8_t op, uint32_t len)
  {
    if (len==0)
      return;
    if (cigar_.empty() || this->op(size()-1)!=op)
      cigar_.push_back(bam_cigar_gen(len, op));
    else
      cigar_[size()-1] += (len<<BAM_CIGAR_SHIFT);
  }


  void Aligner::Cigar::push_front(uint8_t op, uint32_t len)
  {
    if (len==0)
      return;
    if (cigar_.empty() || this->op(0)!=op)
      cigar_.push_front(bam_cigar_gen(len, op));
    else
      cigar_[0] += (len<<BAM_CIGAR_SHIFT);
  }


  uint32_t Aligner::Cigar::query_length(void) const
  {
    return length(1);
  }


  uint32_t Aligner::Cigar::reference_length(void) const
  {
    return length(2);
  }


  size_t Aligner::Cigar::size(void) const
  {
    return cigar_.size();
  }


  uint32_t Aligner::Cigar::operator[](size_t i) const
  {
    assert(i<cigar_.size());
    return cigar_[i];
  }


  std::ostream& operator<<(std::ostream& os, const Aligner::Cigar& cigar)
  {
    for (size_t i=0; i<cigar.size(); ++i)
      os << cigar.oplen(i) << cigar.opchr(i);
    return os;
  }

}}} // of namespace utility, yat, and theplu