source: trunk/yat/utility/SmithWaterman.h @ 3391

#ifndef _theplu_yat_utility_smith_waterman
#define _theplu_yat_utility_smith_waterman

// $Id: SmithWaterman.h 3391 2015-03-16 07:51:09Z peter $

/*
  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 "Aligner.h"
#include "Matrix.h"

#include <boost/concept_check.hpp>

#include <cstddef> // size_t

namespace theplu {
namespace yat {
namespace utility {

  /**
     \since New in yat 0.12
   */
  class SmithWaterman
  {
  public:
    /**
       \brief Constructor

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

    /**
       \brief Constructor

       \param vertical_gap Penalty for extending a vertical gap. A
       vertical gap means consuming a reference element i.e. a
       deletion in query sequence.

       \param open_vertical_gap Open a deletion. Total
       cost for a deletion is open_vertical_gap + N *
       vertical_gap.

       \param horizon_gap Penalty for extending a insertion.

       \param open_horizon_gap Penalty for open an insertion. Total
       penalty for a insertion is open_horizon_gap + N * horizon_gap.
     */
    SmithWaterman(double vertical_gap, double open_vertical_gap,
                  double horizon_gap, double open_horizon_gap);

    /**
       The CIGAR here is slightly different compared to CIGAR in
       Aligner as it uses BAM_CEQUAL and BAM_CDIFF rather than
       BAM_CMATCH, and it also uses BAM_CSOFT_CLIP and ends of CIGAR
       to reflect if whole query was aligned.

       \return CIGAR reflecting latest performed alignment
     */
    const Aligner::Cigar& cigar(void) const;

    /**
       \return position i.e. first element in reference that is
       aligned to the query.
     */
    size_t position(void) const;

    /**
       \return score matrix
     */
    const Matrix& score(void) const;

    /**
       Find best alignment using Smith-Waterman algorithm based on the
       dot-matrix.

       \return smith-waterman score
     */
    double operator()(const Matrix& dot);

    /**
       Calculate a dot-matrix of size query_end-query_begin X
       reference_end-reference_begin where the element ij is 1.0 if
       reference_begin[i] == query_begin[j] and -mismatch otherwise.

       Type Requirements:
       - \c RandomAccessIterator1 is a \readable_iterator
       - \c RandomAccessIterator1 is a \random_access_traversal_iterator
       - \c RandomAccessIterator2 is a \readable_iterator
       - \c RandomAccessIterator2 is a \random_access_traversal_iterator
     */
    template<typename RandomAccessIterator1, typename RandomAccessIterator2>
    double operator()(RandomAccessIterator1 reference_begin,
                      RandomAccessIterator1 reference_end,
                      RandomAccessIterator2 query_begin,
                      RandomAccessIterator2 query_end,
                      double mismatch=0);

  private:
    Aligner aligner_;
    Aligner::Cigar cigar_;
    size_t position_;
    Matrix score_;
  };

  // template implementaion

  template<typename RandomAccessIterator1, typename RandomAccessIterator2>
  double SmithWaterman::operator()(RandomAccessIterator1 reference_begin,
                                   RandomAccessIterator1 reference_end,
                                   RandomAccessIterator2 query_begin,
                                   RandomAccessIterator2 query_end,
                                   double mismatch)
  {
    using boost_concepts::ReadableIterator;
    using boost_concepts::RandomAccessTraversal;
    BOOST_CONCEPT_ASSERT((ReadableIterator<RandomAccessIterator1>));
    BOOST_CONCEPT_ASSERT((RandomAccessTraversal<RandomAccessIterator1>));
    BOOST_CONCEPT_ASSERT((ReadableIterator<RandomAccessIterator2>));
    BOOST_CONCEPT_ASSERT((RandomAccessTraversal<RandomAccessIterator2>));
    Matrix dot(reference_end-reference_begin, query_end-query_begin, -mismatch);
    for (size_t i=0; i<dot.rows(); ++i)
      for (size_t j=0; j<dot.columns(); ++j)
        if (reference_begin[i] == query_begin[j])
          dot(i, j) = 1;
    return this->operator()(dot);
  }

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

#endif