GranthamWeights.cc

/*

  Copyright (C) 2003-2009 Kevin Thornton, krthornt[]@[]uci.edu

  Remove the brackets to email me.

  This file is part of libsequence.

  libsequence 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.

  libsequence 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
  long with libsequence.  If not, see <http://www.gnu.org/licenses/>.

*/

#include <cmath>
#include <cfloat>
#include <Sequence/Grantham.hpp>
#include <Sequence/Translate.hpp>
#include <Sequence/PathwayHelper.hpp>
#include <Sequence/GranthamWeights.hpp>

using std::string;

namespace Sequence
{
  WeightingScheme2::weights2_t GranthamWeights2::operator()(const std::string &codon1, const std::string &codon2,
                                                            GeneticCodes code) const
  {
    Grantham gdist;
    auto intermediates = Intermediates2(codon1,codon2);
    //assign weights to pathways
    //weights are assiged by the total length of each path, as
    //measured by Grantham distances
    double len_path_1 = 0.0, len_path_2 = 0.0;

    string t1 = Translate (codon1.begin(),codon1.end(), code);
    string t2 = Translate (intermediates[0].begin(),
                           intermediates[0].end(), code);
    len_path_1 += gdist ((t1)[0], (t2)[0]);

    t1 = Translate (intermediates[0].begin(),intermediates[0].end(), code);
    t2 = Translate (codon2.begin(),codon2.end(), code);
    len_path_1 += gdist ((t1)[0], (t2)[0]);

    t1 = Translate (codon1.begin(),codon1.end(), code);
    t2 = Translate (intermediates[1].begin(),intermediates[1].end(), code);
    len_path_2 += gdist ((t1)[0], (t2)[0]);

    t1 = Translate (intermediates[1].begin(),intermediates[1].end(), code);
    t2 = Translate (codon2.begin(),codon2.end(), code);
    len_path_2 += gdist ((t1)[0], (t2)[0]);

    //calculate the weights themselves
    //double w_path1 = 0., w_path2 = 0.,
    double w_tot = 0.;
    weights2_t __weights;
    __weights[0]=0.;
    __weights[1]=0.;

    if (fabs(len_path_1-0.) <= DBL_EPSILON && fabs(len_path_2-0.) <= DBL_EPSILON)
      {
        __weights[0] = __weights[1] = 0.5;
        w_tot = 1.;
      }
    else
      {
        __weights[0] = 1. - len_path_1 / (len_path_1 + len_path_2);
        __weights[1] = 1. - len_path_2 / (len_path_1 + len_path_2);
        w_tot = __weights[0] + __weights[1];
      }

    __weights[0] /= w_tot;
    __weights[1] /= w_tot;
    return __weights;
  }

  WeightingScheme3::weights3_t GranthamWeights3::operator()(const std::string &codon1, const std::string &codon2,
                                                            Sequence::GeneticCodes code) const
  {
    Grantham gdist;
    auto intermediates = Intermediates3(codon1,codon2);
    double len_path_1 = 0.0, len_path_2 = 0.0, len_path_3 =
      0., len_path_4 = 0., len_path_5 = 0., len_path_6 = 0.;

    string t1 = Translate (codon1.begin(), codon1.end(),code);
    string t2 = Translate (intermediates[0].begin(), intermediates[0].end(),code);
    double dist = gdist ((t1)[0],(t2)[0]);
    len_path_1 += dist;


    t1 = Translate (intermediates[0].begin(), intermediates[0].end(), code);
    t2 = Translate (intermediates[1].begin(), intermediates[1].end(), code);
    dist = gdist ((t1)[0],(t2)[0]);
    len_path_1 += dist;


    t1 = Translate (intermediates[1].begin(), intermediates[1].end(),code);
    t2 = Translate (codon2.begin(),codon2.end(), code);
    dist = gdist ((t1)[0],(t2)[0]);
    len_path_1 += dist;

    //path2

    t1 = Translate (codon1.begin(),codon1.end(), code);
    t2 = Translate (intermediates[0].begin(),intermediates[0].end(), code);
    dist = gdist ((t1)[0],(t2)[0]);
    len_path_2 += dist;

    t1 = Translate (intermediates[0].begin(),intermediates[0].end(), code);
    t2 = Translate (intermediates[2].begin(),intermediates[2].end(), code);
    dist = gdist ((t1)[0],(t2)[0]);
    len_path_2 += dist;

    t1 = Translate (intermediates[2].begin(),intermediates[2].end(), code);
    t2 = Translate (codon2.begin(),codon2.end(), code);
    dist = gdist ((t1)[0],(t2)[0]);
    len_path_2 += dist;

    //path 3
    t1 = Translate (codon1.begin(),codon1.end(), code);
    t2 = Translate (intermediates[3].begin(),intermediates[3].end(), code);
    dist = gdist ((t1)[0],(t2)[0]);
    len_path_3 += dist;

    t1 = Translate (intermediates[3].begin(), intermediates[3].end(),code);
    t2 = Translate (intermediates[4].begin(), intermediates[4].end(), code);
    dist = gdist ((t1)[0],(t2)[0]);
    len_path_3 += dist;

    t1 = Translate (intermediates[4].begin(), intermediates[4].end(), code);
    t2 = Translate (codon2.begin(),codon2.end(), code);
    dist = gdist ((t1)[0],(t2)[0]);
    len_path_3 += dist;

    //path4
    t1 = Translate (codon1.begin(),codon1.end(), code);
    t2 = Translate (intermediates[3].begin(), intermediates[3].end(), code);
    dist = gdist ((t1)[0],(t2)[0]);
    len_path_4 += dist;

    t1 = Translate (intermediates[3].begin(), intermediates[3].end(), code);
    t2 = Translate (intermediates[5].begin(), intermediates[5].end(), code);
    dist = gdist ((t1)[0],(t2)[0]);
    len_path_4 += dist;

    t1 = Translate (intermediates[5].begin(), intermediates[5].end(), code);
    t2 = Translate (codon2.begin(),codon2.end(), code);
    dist = gdist ((t1)[0],(t2)[0]);
    len_path_4 += dist;

    //path 5
    t1 = Translate (codon1.begin(),codon1.end(), code);
    t2 = Translate (intermediates[6].begin(), intermediates[6].end(), code);
    dist = gdist ((t1)[0],(t2)[0]);
    len_path_5 += dist;

    t1 = Translate (intermediates[6].begin(), intermediates[6].end(), code);
    t2 = Translate (intermediates[7].begin(), intermediates[7].end(), code);
    dist = gdist ((t1)[0],(t2)[0]);
    len_path_5 += dist;

    t1 = Translate (intermediates[7].begin(), intermediates[7].end(), code);
    t2 = Translate (codon2.begin(),codon2.end(), code);
    dist = gdist ((t1)[0],(t2)[0]);
    len_path_5 += dist;

    //path 6
    t1 = Translate (codon1.begin(),codon1.end(), code);
    t2 = Translate (intermediates[6].begin(), intermediates[6].end(), code);
    dist = gdist ((t1)[0],(t2)[0]);
    len_path_6 += dist;

    t1 = Translate (intermediates[6].begin(), intermediates[6].end(), code);
    t2 = Translate (intermediates[8].begin(), intermediates[8].end(), code);
    dist = gdist ((t1)[0],(t2)[0]);
    len_path_6 += dist;

    t1 = Translate (intermediates[8].begin(), intermediates[8].end(), code);
    t2 = Translate (codon2.begin(),codon2.end(), code);
    dist = gdist ((t1)[0],(t2)[0]);
    len_path_6 += dist;

    weights3_t __weights;
    __weights[0] = 1. / len_path_1;
    __weights[1] = 1. /len_path_2;
    __weights[2] = 1. /len_path_3;
    __weights[3] = 1. / len_path_4;
    __weights[4] = 1. /len_path_5;
    __weights[5] = 1./ len_path_6;

    //scale weights to sum to 1
    double w_tot =
      __weights[0] + __weights[1] + __weights[2] + __weights[3] + __weights[4] + __weights[5];
    __weights[0] /= w_tot;
    __weights[1] /= w_tot;
    __weights[2] /= w_tot;
    __weights[3] /= w_tot;
    __weights[4] /= w_tot;
    __weights[5] /= w_tot;
    return __weights;
  }
}