source: trunk/yat/utility/Aligner.h @ 3194

#ifndef _theplu_yat_utility_aligner_
#define _theplu_yat_utility_aligner_

// $Id: Aligner.h 3194 2014-04-22 09:24:56Z peter $

/*
  Copyright (C) 2012 Peter Johansson
  Copyright (C) 2013 Jari Häkkinen, 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 <boost/cstdint.hpp>

#include <cstddef>
#include <iosfwd>
#include <vector>

namespace theplu {
namespace yat {
namespace utility {

  class Matrix;

  /**
     \brief Aligning two sequences

     General class aligning two sequences described in a dot-matrix,
     D, in which \f$ D(i,j) \f$ describes how similar element \c i in
     the first sequence is with element \c j in the second sequence. A
     path through this matrix corresponds to an alignment of the two
     sequences, in which a diagonal step correspoinds to a match
     between corresponding sequence elements, and a vertical or
     horizontal step correspond correspond to inserting a gap into one
     of the sequences.

     \since New in yat 0.9
   */
  class Aligner
  {
  public:
    /**
       Enum used to describe alignment.

       \see alignment
     */
    enum direction { none, right, down, diagonal };

    /**
       \brief Constructor

       Same as Aligner(gap, open_gap, gap, open_gap)
     */
    Aligner(double gap=0, double open_gap=0);

    /**
       \brief constructor

       \param vertical_gap Penalty for extending a vertical gap. A
       vertical gap means alignment consumes an element in second
       element and not, in the first element, i.e., an element has
       been inserted to the second element or equivalently an element
       has been deleted in first sequence.

       \param open_vertical_gap Penalty for open a vertical gap. Total
       cost for a vertical gap is open_vertical_gap + N *
       vertical_gap.

       \param horizon_gap Penalty for extending a horizontal gap

       \param open_horizon_gap Penalty for open a horizontal
       gap. Total penalty for a horizontal gap is open_horizon_gap + N
       * horizon_gap.
     */
    Aligner(double vertical_gap, double open_vertical_gap,
            double horizon_gap, double open_horizon_gap);

    /**
       Initialize a score matrix with \f$ S_{k,0} = -
       \textrm{open\_vertical\_gap} - k * \textrm{vertical\_gap} \f$
       and \f$ S_{0,k} = - \textrm{open\_horizon\_gap} - k *
       \textrm{horizon\_gap} \f$ and align using operator().

       \return most lower right element of score matrix
     */
    double needleman_wunsch(const Matrix& d);

    /**
       Initialize a score matrix with all elements equal to zero,
       align using operator().

       \return max element in score matrix.
     */
    double smith_waterman(const Matrix& d);

    /**
       \brief calculate score matrix based on the \a dot matrix

       The algorithm calculates a maximum \a score matrix as

       \f$ S_{i,j} = \textrm{max} \{ S_{ij}; S_{i,j-1} - F(\textrm{gap}_V);
       S_{i-1,j-1} + D_{i-1,j-1}; S_{i-1,j} - F(\textrm{gap}_H) \} \f$

       where \f$ F(\textrm{gap}) \f$ is gap if there was already a
       gap, and gap + open_gap otherwise.

     */
    void operator()(const Matrix& dot, Matrix& score);

    /**
       If, for example, alignment(i,j) returns Aligner::diagonal,
       implies that score(i,j) was maximized as score(i-1,j-1) +
       dot(i-1, j-1), in other words, element i-1 in first sequence
       was aligned with element j-1 in second sequence.
    */
    const direction& alignment(size_t i, size_t j) const;

    /**
       \since new in yat 0.12
     */
    class Cigar
    {
    public:
      void clear(void);
      void push_back(uint8_t op, uint32_t len=1);
      void push_front(uint8_t op, uint32_t len=1);
      void pop_back(void);
      void pop_front(void);
      uint8_t op(size_t) const;
      char opchr(size_t) const;
      uint32_t oplen(size_t) const;
      uint32_t operator[](size_t) const;
      size_t size(void) const;
    private:
      // using compiler generated copy
      // Cigar(const Cigar& other);
      // Cigar& operator=(const Cigar&);
    };

    /**
       \since new in yat 0.12
     */
    const Cigar cigar(size_t i, size_t j) const;

  private:
    direction& directions(size_t i, size_t j);

    size_t columns_;
    double vertical_gap_;
    double open_vertical_gap_;
    double horizon_gap_;
    double open_horizon_gap_;
    std::vector<direction> alignment_;
  };


  /**
   */
  std::ostream& operator<<(std::ostream& os, const Aligner::Cigar& cigar);

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

#endif