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

#ifndef _theplu_yat_utility_aligner_
#define _theplu_yat_utility_aligner_

// $Id: Aligner.h 3211 2014-05-05 06:38:55Z peter $

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

#include <deque>
#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)
     */
    explicit 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.

       \see SmithWaterman
     */
    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.

       To get the CIGAR to behave as expected, the each row in \a dot
       corresponds to an element in reference sequence and each row in
       \a dot corresponds to an element in query reference. If that is
       transposed, insertions and deletions are swapped in CIGAR.
     */
    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;

    /**
       \brief Compact Idiosyncratic Gapped Alignment Report

       A CIGAR is a representation of an alignment, an array of
       operation/length where the operation and length is packed into
       a uint32_t. So a CIGAR of e.g. '3M1D5M' means alignment is
       three matches, one deletion, and five matches, which is
       represented by an array of length 3.

       <table>
       <tr><td>\c \#define</td><td>value</td><td>char</td><td>type</td>
       <td>consume ref</td><td>consume query</td></tr>
       <tr><td>BAM_CMATCH</td>    <td>0</td><td>M</td><td>3</td>
       <td>\c true</td><td>\c true</td></tr>
       <tr><td>BAM_CINS</td>      <td>1</td><td>I</td><td>1</td>
       <td>\c false</td><td>\c true</td></tr>
       <tr><td>BAM_CDEL</td>      <td>2</td><td>D</td><td>2</td>
       <td>\c true</td><td>\c false</td></tr>
       <tr><td>BAM_CREF_SKIP</td> <td>3</td><td>N</td><td>2</td>
       <td>\c true</td><td>\c false</td></tr>
       <tr><td>BAM_CSOFT_CLIP</td><td>4</td><td>S</td><td>1</td>
       <td>\c false</td><td>\c true</td></tr>
       <tr><td>BAM_CHARD_CLIP</td><td>5</td><td>H</td><td>0</td>
       <td>\c false</td><td>\c false</td></tr>
       <tr><td>BAM_CPAD</td>      <td>6</td><td>P</td><td>0</td>
       <td>\c false</td><td>\c false</td></tr>
       <tr><td>BAM_CEQUAL</td>    <td>7</td><td>=</td><td>3</td>
       <td>\c true</td><td>\c true</td></tr>
       <tr><td>BAM_CDIFF</td>     <td>8</td><td>X</td><td>3</td>
       <td>\c true</td><td>\c true</td></tr>
       <tr><td>BAM_CBACK</td>     <td>9</td><td>B</td><td>0</td>
       <td>\c false</td><td>\c false</td></tr>
       </table>

       \since new in yat 0.12
     */
    class Cigar
    {
    public:
      /**
         Erases all elements.
       */
      void clear(void);

      /**
         See table in class documntation.

         \return operation part of element \a i
       */
      uint8_t op(size_t i) const;

      /**
         Translate the op(i) to a descriptive character of one of the
         following "MIDNSHP=XB". See table in class documentation.

         \return character describing operation i
       */
      char opchr(size_t i) const;

      /**
         \return length of element i
       */
      uint32_t oplen(size_t i) const;

      /**
         Erase last \a len operation.
       */
      void pop_back(uint32_t len=1);

      /**
         Erase first \a len operation.
       */
      void pop_front(uint32_t len=1);

      /**
         Add an operation \a op of length \a len at back of cigar. If
         last operation equals \a op, length of that element is
         changed and the size() is not modified.
       */
      void push_back(uint8_t op, uint32_t len=1);

      /**
         Add an operation \a op of length \a len at front of cigar.
       */
      void push_front(uint8_t op, uint32_t len=1);

      /**
         \return cigar element \a i
       */
      uint32_t operator[](size_t i) const;

      /**
         \return number of elements in cigar
       */
      size_t size(void) const;
    private:
      std::deque<uint32_t> cigar_;

      // 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_;
  };


  /**
     Output e.g. '5S44M5I98M' if the cigar has four elements. Each
     element is output as the length followed by the character
     descibing the operation (see opchr).

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

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

#endif