//------------------------------------------------------------------------------ // This file is part of the OpenStructure project <www.openstructure.org> // // Copyright (C) 2008-2011 by the OpenStructure authors // // This library is free software; you can redistribute it and/or modify it under // the terms of the GNU Lesser General Public License as published by the Free // Software Foundation; either version 3.0 of the License, or (at your option) // any later version. // This 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 Lesser General Public License for more // details. // // You should have received a copy of the GNU Lesser General Public License // along with this library; if not, write to the Free Software Foundation, Inc., // 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA //------------------------------------------------------------------------------ #include <ost/log.hh> #include "impl/align_impl.hh" #include "global_align.hh" namespace ost { namespace seq { namespace alg { void Traceback(impl::AlnMat& mat, int max_i, int max_j, const ConstSequenceHandle& s1, const ConstSequenceHandle& s2, AlignmentList& alignments) { int i = (max_i > 0 ? max_i-1 : 0); int j = (max_j > 0 ? max_j-1 : 0); String aln_str1; String aln_str2; while (i > 0 || j > 0) { impl::SetRoute(mat, i, j, s1, s2, aln_str1, aln_str2); } impl::StoreStrAsAln(aln_str1, aln_str2, s1, s2, i, j, alignments); } AlignmentList GlobalAlign(const ConstSequenceHandle& s1, const ConstSequenceHandle& s2, alg::SubstWeightMatrixPtr& subst, int gap_open,int gap_ext) { impl::AlnMat mat(s1.GetLength()+1, s2.GetLength()+1); // init first column of matrix if (mat.GetHeight() > 1) { mat(0, 1).score = mat(0, 0).score + gap_open; mat(0, 1).from = impl::INS1; } for (int j = 2; j < mat.GetHeight(); ++j) { mat(0, j).score = mat(0, j-1).score + gap_ext; mat(0, j).from = impl::INS1; } // fill the alignment matrix for (int i = 0; i < mat.GetWidth()-1; ++i) { mat(i+1, 0).score = mat(i, 0).score + (mat(i, 0).from==impl::INS2 ? gap_ext : gap_open); mat(i+1,0).from = impl::INS2; for (int j = 0; j < mat.GetHeight()-1; ++j) { char c1=s1[i]; char c2=s2[j]; short weight=subst->GetWeight(c1, c2); int diag=weight+mat(i, j).score; int ins1=mat(i+1, j).score +(mat(i+1, j).from==impl::INS1 ? gap_ext : gap_open); int ins2=mat(i, j+1).score +(mat(i, j+1).from==impl::INS2 ? gap_ext : gap_open); if (diag>=ins1) { if (diag>=ins2) { mat(i+1, j+1).score=diag; mat(i+1, j+1).from=impl::DIAG; } else { mat(i+1, j+1).score=ins2; mat(i+1, j+1).from=impl::INS2; } } else if (ins1>ins2) { mat(i+1, j+1).score=ins1; mat(i+1, j+1).from=impl::INS1; } else { mat(i+1, j+1).score=ins2; mat(i+1, j+1).from=impl::INS2; } } } // write traceback matrix in debug mode #if !defined(NDEBUG) DbgWriteAlnMatrix(mat,s1, s2); #endif AlignmentList alignments; Traceback(mat, mat.GetWidth(), mat.GetHeight(), s1, s2, alignments); LOG_DEBUG(alignments.back().ToString(80)); return alignments; } }}}