Added functions to predict contacts from MSAs.

modules/seq/alg/pymod/__init__.py

@@ -168,3 +168,37 @@ def AlignmentFromChainView(chain, handle_seq_name='handle',
       s1.Append('-')
     idx0+=1
   return seq.CreateAlignment(s0, s1)
+
+def PredictContacts(ali):
+  """
+  Predicts contacts from a multiple sequence alignment using a combination
+  of Mutual Information (*MI*) and the Contact Substitution Score (*CoEvoSc*).
+  MI is calculated with the APC and small number corrections as well as with a 
+  transformation into Z-scores. The *CoEvoSc* is calculated using the default 
+  PairSubstWeightMatrix (see seq.alg.LoadDefaultPairSubstWeightMatrix).
+  The final score for a pair of columns *(i,j)* of **ali** is obtained from:
+  
+  Sc(i,j)=MI(i,j)exp(CoEvoSc(i,j))      if *(i,j)* >=0
+  
+  Sc(i,j)=MI(i,j)exp(1-CoEvoSc(i,j))    if *(i,j)* <0
+
+  :param ali: The multiple sequence alignment
+  :type ali: :class:`~ost.seq.AlignmentHandle`
+  """
+  import math
+  from ost import seq
+  if not type(ali)==type(seq.AlignmentHandle()):
+    print "Parameter should be an AlignmentHandle"
+    return
+  mi=CalculateMutualInformation(ali)
+  CoEvoSc=CalculateContactSubstitutionScore(ali)
+  ncol=ali.GetLength()
+  for i in range(ncol):
+    for j in range(ncol):
+      if mi.matrix[i][j]>=0:
+        mi.matrix[i][j]=mi.matrix[i][j]*(math.exp(CoEvoSc.matrix[i][j]))
+      else:
+        mi.matrix[i][j]=mi.matrix[i][j]*(math.exp(1.-CoEvoSc.matrix[i][j]))
+  mi.RefreshSortedIndices()
+  return mi
+

modules/seq/alg/pymod/wrap_seq_alg.cc

@@ -29,6 +29,10 @@
 #include <ost/seq/alg/global_align.hh>
 #include <ost/seq/alg/semiglobal_align.hh>
 #include <ost/seq/alg/entropy.hh>
+#include <ost/seq/alg/pair_subst_weight_matrix.hh>
+#include <ost/seq/alg/contact_weight_matrix.hh>
+#include <ost/seq/alg/contact_prediction_score.hh>
+
 using namespace boost::python;
 using namespace ost::seq;
 using namespace ost::seq::alg;
@@ -68,6 +72,32 @@ BOOST_PYTHON_MODULE(_ost_seq_alg)
       arg("gap_open")=-5, arg("gap_ext")=-2));
   def("ShannonEntropy", &ShannonEntropy, (arg("aln"), arg("ignore_gaps")=true));
 
+  class_<PairSubstWeightMatrix>("PairSubstWeightMatrix",init< std::vector <std::vector <std::vector <std::vector <Real> > > >,std::vector <char> >())
+    .def(init<>())
+    .def_readonly("weights",&PairSubstWeightMatrix::weights)
+    .def_readonly("aa_list",&PairSubstWeightMatrix::aa_list)
+  ;
+
+  class_<ContactWeightMatrix>("ContactWeightMatrix",init< std::vector <std::vector <Real> >,std::vector <char> >())
+    .def(init<>())
+    .def_readonly("weights",&ContactWeightMatrix::weights)
+    .def_readonly("aa_list",&ContactWeightMatrix::aa_list)
+  ;
+
+  class_<ContactPredictionScoreResult>("ContactPredictionScoreResult",no_init)
+    .def_readonly("matrix",&ContactPredictionScoreResult::matrix)
+    .def_readonly("sorted_indices",&ContactPredictionScoreResult::sorted_indices)
+    .def("RefreshSortedIndices",&ContactPredictionScoreResult::RefreshSortedIndices)
+    .def("GetScore",&ContactPredictionScoreResult::GetScore,(arg("i"),arg("j")))
+    .def("SetScore",&ContactPredictionScoreResult::SetScore,(arg("i"),arg("j"),arg("score")))
+  ;
+  def("CalculateMutualInformation", &CalculateMutualInformation,(arg("aln"),arg("w")=LoadConstantContactWeightMatrix(),
+                                                        arg("apc_correction")=true,arg("zpx_transformation")=true,arg("small_number_correction")=0.05));
+  def("CalculateContactScore", &CalculateContactScore,(arg("aln"), arg("w")=LoadDefaultContactWeightMatrix()));
+  def("CalculateContactSubstitutionScore", &CalculateContactSubstitutionScore,(arg("aln"), arg("ref_seq_index")=0, arg("w")=LoadDefaultPairSubstWeightMatrix()));
+  def("LoadDefaultContactWeightMatrix",LoadDefaultContactWeightMatrix);
+  def("LoadConstantContactWeightMatrix",LoadConstantContactWeightMatrix);
+  def("LoadDefaultPairSubstWeightMatrix",LoadDefaultPairSubstWeightMatrix);
 
   scope mat_scope = class_<SubstWeightMatrix, SubstWeightMatrixPtr>("SubstWeightMatrix", init<>())
                       .def("GetWeight", &SubstWeightMatrix::GetWeight)

modules/seq/alg/src/CMakeLists.txt

@@ -11,6 +11,11 @@ conservation.hh
 local_align.hh
 global_align.hh
 semiglobal_align.hh
+pair_subst_weight_matrix.hh
+contact_weight_matrix.hh
+contact_prediction_score.hh
+data/default_pair_subst_weight_matrix.hh
+data/default_contact_weight_matrix.hh
 )
 
 set(OST_SEQ_ALG_SOURCES
@@ -23,7 +28,10 @@ sequence_identity.cc
 sequence_similarity.cc
 ins_del.cc
 conservation.cc
+contact_weight_matrix.cc
 subst_weight_matrix.cc
+pair_subst_weight_matrix.cc
+contact_prediction_score.cc
 )
 
 module(NAME seq_alg HEADER_OUTPUT_DIR ost/seq/alg SOURCES ${OST_SEQ_ALG_SOURCES}

modules/seq/alg/src/contact_prediction_score.cc

+//------------------------------------------------------------------------------
+// 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/seq/aligned_column.hh>
+#include <ost/seq/alignment_handle.hh>
+#include <ost/seq/alg/contact_prediction_score.hh>
+#include <math.h>
+
+namespace {
+  typedef std::pair<float, std::pair<int,int> > val_indices_pair;
+
+  bool pair_comparator ( const val_indices_pair& l, const val_indices_pair& r)
+    { return l.first > r.first; }
+}
+
+namespace ost { namespace seq { namespace alg {
+
+void ContactPredictionScoreResult::RefreshSortedIndices(){
+  int n_col=this->matrix.size();
+  std::vector <val_indices_pair> sorting_vec;
+  for (int col1=0; col1<n_col; ++col1) {
+      for (int col2=col1+1; col2<n_col; ++col2) {
+        std::pair <int,int> indices (col1,col2);
+        val_indices_pair vi (this->matrix[col1][col2],indices);
+        sorting_vec.push_back(vi);
+      }
+    }
+  std::sort(sorting_vec.begin(),sorting_vec.end(),pair_comparator);
+  int n_indices=sorting_vec.size();
+  std::vector <std::pair <int,int> > sorted_indices;
+  for (int i=0; i<n_indices; ++i) {
+    sorted_indices.push_back(sorting_vec[i].second);
+  }
+  this->sorted_indices=sorted_indices;
+}
+
+ContactPredictionScoreResult::ContactPredictionScoreResult(std::vector <std::vector <Real> > scores){
+  this->matrix=scores;
+  this->RefreshSortedIndices();
+}
+
+Real ContactPredictionScoreResult::GetScore(int i, int j){
+    return this->matrix[i][j];
+  }
+
+void ContactPredictionScoreResult::SetScore(int i, int j, Real score){
+    this->matrix[i][j]=score;
+  }
+
+
+ContactPredictionScoreResult CalculateContactSubstitutionScore(const AlignmentHandle& aln, int ref_seq_index, PairSubstWeightMatrix w){
+  int n_col=aln.GetLength();
+  int n_row=aln.GetCount();
+  //We transform the alignment in a list of ints
+  std::vector< std::vector <int> > ali(n_col, std::vector <int> (n_row) );
+  int naa=w.aa_list.size();
+  for (int i=0; i<n_col; ++i) {
+    for (int j=0; j<n_row; ++j) {
+      std::map<char,int>::iterator iter=w.aa_dict.find(aln[i][j]);
+      if (iter!=w.aa_dict.end()){
+        ali[i][j]=w.aa_dict[aln[i][j]];
+      }
+      else ali[i][j]=naa;
+    }
+  }
+  std::vector< std::vector <Real> > contact_subst_score(n_col, std::vector <Real> (n_col,0) );
+  for (int col1=0; col1<n_col; ++col1) {
+    std::vector<int> c1=ali[col1];
+    int i1=c1[ref_seq_index];
+    if (i1==naa){
+      for (int col2=col1+1; col2<n_col; ++col2) {
+        contact_subst_score[col1][col2]=0.0;
+        contact_subst_score[col2][col1]=0.0;
+      }
+      continue;
+    }
+    for (int col2=col1+1; col2<n_col; ++col2) {
+      std::vector<int> c2=ali[col2];
+      int i2=c2[ref_seq_index];
+      if (i2==naa){
+        contact_subst_score[col1][col2]=0.0;
+        contact_subst_score[col2][col1]=0.0;
+        continue;
+        }
+      Real score=0.0;
+      for (int i=0 ; i<n_row; ++i){
+        if (i==ref_seq_index) continue;
+        int i3=c1[i], i4=c2[i];
+        if (i3==naa || i4==naa){
+          score+=0.0;
+          continue;
+        }
+        score+=w.weights[i1][i2][i3][i4];
+      }
+    contact_subst_score[col1][col2]=score/(n_row-1.);
+    contact_subst_score[col2][col1]=score/(n_row-1.);
+    }
+  }
+  return ContactPredictionScoreResult(contact_subst_score);
+}
+
+ContactPredictionScoreResult CalculateContactScore(const AlignmentHandle& aln, ContactWeightMatrix weights){
+  int n_col=aln.GetLength();
+  int n_row=aln.GetCount();
+  //We transform the alignment in a list of ints
+  std::vector< std::vector <int> > ali(n_col, std::vector <int> (n_row) );
+  int naa=weights.aa_list.size();
+  for (int i=0; i<n_col; ++i) {
+    for (int j=0; j<n_row; ++j) {
+      std::map<char,int>::iterator iter=weights.aa_dict.find(aln[i][j]);
+      if (iter!=weights.aa_dict.end()){
+        ali[i][j]=weights.aa_dict[aln[i][j]];
+      }
+      else ali[i][j]=naa;
+    }
+  }
+
+  std::vector< std::vector <Real> > mi(n_col, std::vector <Real> (n_col,0) );
+
+  for (int col1=0; col1<n_col; ++col1) {
+    std::vector<int> c1=ali[col1];
+    for (int col2=col1+1; col2<n_col; ++col2) {
+      float score=0.0;
+      for (int i=0 ; i<n_row; ++i){
+        int aa1=c1[i];
+        int aa2=ali[col2][i];
+        if (aa1==naa || aa2==naa) continue;
+        score+=weights.weights[aa1][aa2];
+      }
+      mi[col1][col2]=score/n_row;
+      mi[col2][col1]=score/n_row;
+    }
+  }
+  return ContactPredictionScoreResult(mi);
+}
+
+
+ContactPredictionScoreResult CalculateMutualInformation(const AlignmentHandle& aln, ContactWeightMatrix weights,
+                                  bool apc_correction,bool zpx_transformation, float small_number_correction){
+  int n_col=aln.GetLength();
+  int n_row=aln.GetCount();
+  //We transform the alignment in a list of ints
+  std::vector< std::vector <int> > ali(n_col, std::vector <int> (n_row) );
+  int naa=weights.aa_list.size();
+  for (int i=0; i<n_col; ++i) {
+    for (int j=0; j<n_row; ++j) {
+      std::map<char,int>::iterator iter=weights.aa_dict.find(aln[i][j]);
+      if (iter!=weights.aa_dict.end()){
+        ali[i][j]=weights.aa_dict[aln[i][j]];
+      }
+      else ali[i][j]=naa;
+    }
+  }
+
+  std::vector< std::vector <Real> > mi(n_col, std::vector <Real> (n_col,0) );
+
+  for (int col1=0; col1<n_col; ++col1) {
+    std::vector<int> c1=ali[col1];
+    for (int col2=col1+1; col2<n_col; ++col2) {
+      std::vector<float> p1(naa,naa*small_number_correction),p2(naa,naa*small_number_correction);
+      std::vector<std::vector<float> > p12(naa, std::vector<float> (naa,small_number_correction));
+      for (int i=0 ; i<n_row; ++i){
+        int aa1=c1[i];
+        int aa2=ali[col2][i];
+        if (aa1==naa || aa2==naa) continue;
+        p1[aa1]++;
+        p2[aa2]++;
+        p12[aa1][aa2]++;
+      }
+      float n=naa*small_number_correction;
+      for (int aa1=0; aa1!=naa; ++aa1){
+        n+=p1[aa1];
+      }
+      float score=0.0;
+      for (int aa1=0; aa1!=naa; ++aa1){
+        if(p1[aa1]==0.0) continue;
+        for (int aa2=0; aa2!=naa; ++aa2){
+          if(p12[aa1][aa2]==0.0 || p2[aa2]==0.0) continue;
+          score+=weights.weights[aa1][aa2]*p12[aa1][aa2]/n*std::log(n*p12[aa1][aa2]/(p1[aa1]*p2[aa2]));
+        }
+      }
+      mi[col1][col2]=score;
+      mi[col2][col1]=score;
+    }
+  }
+  if (apc_correction==true){
+    float average_mi=0.;
+    std::vector <float> mi_col(n_col,0.);
+    for (int col1=0; col1<n_col; ++col1) {
+      for (int col2=0; col2<n_col; ++col2) {
+        average_mi+=mi[col1][col2];
+        mi_col[col1]+=mi[col1][col2];
+      }
+    }
+    for (int col1=0; col1<n_col; ++col1) {
+      for (int col2=col1+1; col2<n_col; ++col2) {
+      mi[col1][col2]-=mi_col[col1]*mi_col[col2]/average_mi;
+      mi[col2][col1]-=mi_col[col2]*mi_col[col1]/average_mi;
+      }
+    }
+  }
+  if (zpx_transformation==true){
+    std::vector <float> mi_col(n_col,0.);
+    for (int col1=0; col1<n_col; ++col1) {
+      for (int col2=0; col2<n_col; ++col2) {
+        mi_col[col1]+=mi[col1][col2];
+      }
+    }
+    for (int col1=0; col1<n_col; ++col1) {
+        mi_col[col1]/=n_col;
+    }
+    std::vector <float> mi_std_col(n_col,0.);
+    for (int col1=0; col1<n_col; ++col1) {
+      for (int col2=0; col2<n_col; ++col2) {
+        mi_std_col[col1]+=(mi[col1][col2]-mi_col[col1])*(mi[col1][col2]-mi_col[col1]);
+      }
+    }
+    for (int col1=0; col1<n_col; ++col1) {
+        mi_std_col[col1]=sqrt(mi_std_col[col1]/n_col);
+    }
+    for (int col1=0; col1<n_col; ++col1) {
+      for (int col2=col1+1; col2<n_col; ++col2) {
+      mi[col1][col2]=(mi[col1][col2]-mi_col[col1])/mi_std_col[col1]+(mi[col1][col2]-mi_col[col2])/mi_std_col[col2];
+      mi[col2][col1]=mi[col1][col2];
+      }
+    }
+  }  
+  return ContactPredictionScoreResult(mi);
+}
+
+}}}
+

modules/seq/alg/src/contact_prediction_score.hh

+//------------------------------------------------------------------------------
+// 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
+//------------------------------------------------------------------------------
+#ifndef OST_SEQ_ALG_CONTACT_PREDICTION_SCORE_HH
+#define OST_SEQ_ALG_CONTACT_PREDICTION_SCORE_HH
+#include <ost/seq/alignment_handle.hh>
+#include <ost/seq/alg/module_config.hh>
+#include <ost/seq/alg/pair_subst_weight_matrix.hh>
+#include <ost/seq/alg/contact_weight_matrix.hh>
+
+/*
+  Authors: Niklaus Johner
+  Functions and objects used to calculate different scores for predicting contacts
+  from a multiple sequence alignment (MSA). Several scores are available:
+  1. Mutual information. The average product correction (apc_correction) 
+  adds a correction term to account for common phylogenetic origin and inherent
+  entropy of each column. There is also a correction for small number of sequences 
+  in the alignment and the possibility to use weighted mutual information.
+  2. ContactScore. For each residue pair, this score calculates the average interaction
+  energy for the amino-acid pairs in the corresponding columns of the MSA
+  3. ContactSubstitutionScore. For each residue pair, this score calculates 
+  the average energy of the substitutions of amino-acid pairs in the corresponding
+  columns of the MSA.
+ */
+namespace ost { namespace seq { namespace alg {
+
+struct DLLEXPORT_OST_SEQ_ALG ContactPredictionScoreResult {
+  std::vector <std::vector <Real> > matrix;
+  std::vector <std::pair <int,int> > sorted_indices;
+  void SetScore(int i, int j, Real score);
+  Real GetScore(int i, int j);
+  void RefreshSortedIndices();
+  ContactPredictionScoreResult(std::vector <std::vector <Real> >);
+};
+
+ContactPredictionScoreResult DLLEXPORT_OST_SEQ_ALG CalculateContactScore(const AlignmentHandle& aln, ContactWeightMatrix w=LoadDefaultContactWeightMatrix());
+ContactPredictionScoreResult DLLEXPORT_OST_SEQ_ALG CalculateContactSubstitutionScore(const AlignmentHandle& aln, int ref_seq_index=0, PairSubstWeightMatrix w=LoadDefaultPairSubstWeightMatrix());
+ContactPredictionScoreResult DLLEXPORT_OST_SEQ_ALG CalculateMutualInformation(const AlignmentHandle& aln, ContactWeightMatrix w=LoadConstantContactWeightMatrix(),bool apc_correction=true, 
+                            bool zpx_transformation=true,float small_number_correction=0.05);
+
+}}}
+#endif

modules/seq/alg/src/contact_weight_matrix.cc

+//------------------------------------------------------------------------------
+// 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/info/info.hh>
+#include <ost/seq/alg/contact_weight_matrix.hh>
+#include <ost/seq/alg/default_contact_weight_matrix.hh>
+
+namespace {
+std::map<char,int> init_dict(std::vector <char> aal){
+  std::map<char,int> aa_dict;
+  for (int i=0;i!=aal.size();++i){
+    aa_dict[aal[i]]=i;
+  }
+  return aa_dict;
+}
+}
+
+namespace ost { namespace seq { namespace alg {
+
+ContactWeightMatrix LoadConstantContactWeightMatrix(){
+  std::vector <std::vector<Real> > w(21, std::vector<Real>(21,1.0));
+  std::vector <char> aal(21);
+  for (int i=0;i!=21;i++){
+    aal[i]=RAW_CONTACT_WEIGHT_MATRIX_RES_LIST[i];
+  }
+  for (int i1=0;i1!=21;i1++){
+    w[20][i1]=0.0;
+    w[i1][20]=0.0;
+  }
+  return ContactWeightMatrix(w,aal);
+}
+
+ContactWeightMatrix LoadDefaultContactWeightMatrix(){
+  std::vector <std::vector<Real> > w(21, std::vector<Real>(21));
+  std::vector <char> aal(21);
+  for (int i=0;i!=21;i++){
+    aal[i]=RAW_CONTACT_WEIGHT_MATRIX_RES_LIST[i];
+  }
+  for (int i1=0;i1!=21;i1++){
+    for (int i2=0;i2!=21;i2++){
+      w[i1][i2]=RAW_CONTACT_WEIGHT_MATRIX[i1][i2];
+    }
+  }
+  return ContactWeightMatrix(w,aal);
+}
+
+ContactWeightMatrix::ContactWeightMatrix(){
+  std::vector <std::vector <Real> > w;
+  std::vector <char> aal;
+  this->aa_dict=init_dict(this->aa_list);
+  this->weights=w;
+}
+
+ContactWeightMatrix::ContactWeightMatrix(std::vector <std::vector <Real> > w,std::vector <char> aal){
+  this->aa_list=aal;
+  this->naa_=this->aa_list.size();
+  this->aa_dict=init_dict(this->aa_list);
+  this->weights=w;
+}
+
+
+}}}
+

modules/seq/alg/src/contact_weight_matrix.hh

+//------------------------------------------------------------------------------
+// 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
+//------------------------------------------------------------------------------
+#ifndef OST_SEQ_CONTACT_WEIGHT_MATRIX_HH
+#define OST_SEQ_CONTACT_WEIGHT_MATRIX_HH
+
+#include <map>
+#include <ost/seq/alg/module_config.hh>
+
+/*
+  Author: Niklaus Johner
+ */
+namespace ost { namespace seq { namespace alg {
+
+  
+/// \brief matrix to weight pairs of amino-acids.
+struct DLLEXPORT_OST_SEQ_ALG ContactWeightMatrix {
+  std::vector <std::vector <Real> > weights;
+  std::vector <char> aa_list;
+  std::map<char,int> aa_dict;
+  int naa_;
+
+  ContactWeightMatrix();
+  ContactWeightMatrix(std::vector <std::vector <Real> >,std::vector <char>);
+};
+
+/// \brief statistical potential matrix containing interaction pseudo energies
+ContactWeightMatrix DLLEXPORT_OST_SEQ_ALG LoadDefaultContactWeightMatrix();
+/// \brief weight of 1 for all amino-acid pairs and 0 for gaps.
+ContactWeightMatrix DLLEXPORT_OST_SEQ_ALG LoadConstantContactWeightMatrix();
+
+}}}
+
+#endif
\ No newline at end of file

modules/seq/alg/src/data/default_contact_weight_matrix.hh

+namespace ost { namespace seq { namespace alg {
+
+char RAW_CONTACT_WEIGHT_MATRIX_RES_LIST[21]={'D', 'E', 'R', 'K', 'H', 'S', 'T', 'N', 'Q', 'G', 'P', 'Y', 'W', 'V', 'I', 'L', 'M', 'F', 'A', 'C','-'};
+
+Real RAW_CONTACT_WEIGHT_MATRIX[21][21]=
+    {
+      { 0.0168, 0.0000, 0.4730, 0.3602, 0.3627, 0.1650, 0.1716, 0.2177, 0.2272, 0.0770, 0.0733, 0.3501, 0.3262, 0.1228, 0.1393, 0.1263, 0.1792, 0.2105, 0.0757, 0.1801, 0.0000},
+      { 0.0000, 0.0676, 0.4996, 0.3966, 0.3357, 0.1594, 0.1866, 0.1958, 0.2390, 0.0333, 0.0945, 0.3699, 0.3719, 0.1753, 0.1948, 0.1866, 0.2057, 0.2506, 0.1179, 0.1688, 0.0000},
+      { 0.4730, 0.4996, 0.3499, 0.1961, 0.3666, 0.2970, 0.3204, 0.3266, 0.3599, 0.2218, 0.2777, 0.5060, 0.5262, 0.3196, 0.3532, 0.3577, 0.3560, 0.4251, 0.2489, 0.3517, 0.0000},
+      { 0.3602, 0.3966, 0.1961, 0.1128, 0.2100, 0.1731, 0.2061, 0.2093, 0.2493, 0.0799, 0.0802, 0.3928, 0.3725, 0.2150, 0.2331, 0.2226, 0.2164, 0.2845, 0.1383, 0.2104, 0.0000},
+      { 0.3627, 0.3357, 0.3666, 0.2100, 0.5496, 0.3269, 0.3494, 0.3122, 0.3180, 0.2106, 0.2629, 0.5362, 0.5588, 0.3475, 0.3754, 0.3618, 0.4193, 0.4679, 0.2361, 0.4789, 0.0000},
+      { 0.1650, 0.1594, 0.2970, 0.1731, 0.3269, 0.1866, 0.2143, 0.2210, 0.2315, 0.1203, 0.1250, 0.3501, 0.3960, 0.2381, 0.2586, 0.2358, 0.2758, 0.3335, 0.1628, 0.3010, 0.0000},
+      { 0.1716, 0.1866, 0.3204, 0.2061, 0.3494, 0.2143, 0.2894, 0.2553, 0.2783, 0.1418, 0.1529, 0.3962, 0.4310, 0.3452, 0.3660, 0.3237, 0.3417, 0.3951, 0.2073, 0.3458, 0.0000},
+      { 0.2177, 0.1958, 0.3266, 0.2093, 0.3122, 0.2210, 0.2553, 0.3119, 0.2868, 0.1666, 0.1810, 0.3706, 0.3994, 0.2195, 0.2391, 0.2168, 0.2834, 0.3160, 0.1618, 0.2983, 0.0000},
+      { 0.2272, 0.2390, 0.3599, 0.2493, 0.3180, 0.2315, 0.2783, 0.2868, 0.3159, 0.1473, 0.2187, 0.4208, 0.4769, 0.2727, 0.3014, 0.2989, 0.3234, 0.3733, 0.2010, 0.2982, 0.0000},
+      { 0.0770, 0.0333, 0.2218, 0.0799, 0.2106, 0.1203, 0.1418, 0.1666, 0.1473, 0.0990, 0.0727, 0.3054, 0.3354, 0.1831, 0.1984, 0.1959, 0.2488, 0.2887, 0.1321, 0.2889, 0.0000},
+      { 0.0733, 0.0945, 0.2777, 0.0802, 0.2629, 0.1250, 0.1529, 0.1810, 0.2187, 0.0727, 0.1358, 0.4454, 0.4871, 0.1879, 0.2176, 0.2163, 0.2738, 0.3492, 0.1143, 0.2611, 0.0000},
+      { 0.3501, 0.3699, 0.5060, 0.3928, 0.5362, 0.3501, 0.3962, 0.3706, 0.4208, 0.3054, 0.4454, 0.6013, 0.6686, 0.5248, 0.5702, 0.5625, 0.5797, 0.6506, 0.3876, 0.5335, 0.0000},
+      { 0.3262, 0.3719, 0.5262, 0.3725, 0.5588, 0.3960, 0.4310, 0.3994, 0.4769, 0.3354, 0.4871, 0.6686, 0.7438, 0.5554, 0.6133, 0.6157, 0.6450, 0.7208, 0.4170, 0.6399, 0.0000},
+      { 0.1228, 0.1753, 0.3196, 0.2150, 0.3475, 0.2381, 0.3452, 0.2195, 0.2727, 0.1831, 0.1879, 0.5248, 0.5554, 0.5436, 0.5726, 0.5340, 0.4934, 0.5834, 0.3537, 0.4786, 0.0000},
+      { 0.1393, 0.1948, 0.3532, 0.2331, 0.3754, 0.2586, 0.3660, 0.2391, 0.3014, 0.1984, 0.2176, 0.5702, 0.6133, 0.5726, 0.6365, 0.5952, 0.5478, 0.6388, 0.3921, 0.5096, 0.0000},
+      { 0.1263, 0.1866, 0.3577, 0.2226, 0.3618, 0.2358, 0.3237, 0.2168, 0.2989, 0.1959, 0.2163, 0.5625, 0.6157, 0.5340, 0.5952, 0.5835, 0.5272, 0.6279, 0.3691, 0.4760, 0.0000},
+      { 0.1792, 0.2057, 0.3560, 0.2164, 0.4193, 0.2758, 0.3417, 0.2834, 0.3234, 0.2488, 0.2738, 0.5797, 0.6450, 0.4934, 0.5478, 0.5272, 0.5920, 0.6433, 0.3527, 0.5070, 0.0000},
+      { 0.2105, 0.2506, 0.4251, 0.2845, 0.4679, 0.3335, 0.3951, 0.3160, 0.3733, 0.2887, 0.3492, 0.6506, 0.7208, 0.5834, 0.6388, 0.6279, 0.6433, 0.7323, 0.4201, 0.5920, 0.0000},
+      { 0.0757, 0.1179, 0.2489, 0.1383, 0.2361, 0.1628, 0.2073, 0.1618, 0.2010, 0.1321, 0.1143, 0.3876, 0.4170, 0.3537, 0.3921, 0.3691, 0.3527, 0.4201, 0.2846, 0.3407, 0.0000},
+      { 0.1801, 0.1688, 0.3517, 0.2104, 0.4789, 0.3010, 0.3458, 0.2983, 0.2982, 0.2889, 0.2611, 0.5335, 0.6399, 0.4786, 0.5096, 0.4760, 0.5070, 0.5920, 0.3407, 1.0000, 0.0000},
+      { 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000}};
+
+}}}
\ No newline at end of file

modules/seq/alg/src/pair_subst_weight_matrix.cc

+//------------------------------------------------------------------------------
+// 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/info/info.hh>
+#include <ost/seq/alg/pair_subst_weight_matrix.hh>
+#include <ost/seq/alg/default_pair_subst_weight_matrix.hh>
+#include <algorithm>
+
+namespace {
+
+std::map<char,int> init_dict(std::vector <char> aal){
+  std::map<char,int> aa_dict;
+  for (int i=0;i!=aal.size();++i){
+    aa_dict[aal[i]]=i;
+  }
+  return aa_dict;
+}
+}
+
+namespace ost { namespace seq { namespace alg {
+
+PairSubstWeightMatrix LoadDefaultPairSubstWeightMatrix(){
+  std::vector <std::vector<std::vector<std::vector<Real> > > > w(20, std::vector<std::vector<std::vector<Real> > >(20,std::vector<std::vector<Real> >(20,std::vector<Real>(20))));
+  std::vector <char> aal(20);
+  for (int i=0;i!=20;i++){
+    aal[i]=RAW_PAIR_SUBST_WEIGHT_MATRIX_RES_LIST[i];
+  }
+  for (int i1=0;i1!=20;i1++){
+    for (int i2=0;i2!=20;i2++){
+      for (int i3=0;i3!=20;i3++){
+        for (int i4=0;i4!=20;i4++){
+          w[i1][i2][i3][i4]=RAW_PAIR_SUBST_WEIGHT_MATRIX[i1][i2][i3][i4];
+        }
+      }
+    }
+  }
+  return PairSubstWeightMatrix(w,aal);
+}
+
+
+PairSubstWeightMatrix::PairSubstWeightMatrix(){
+  std::vector <std::vector <std::vector <std::vector <Real> > > > w;
+  std::vector <char> aal;
+  this->aa_dict=init_dict(this->aa_list);
+  this->weights=w;
+}
+
+PairSubstWeightMatrix::PairSubstWeightMatrix(std::vector <std::vector <std::vector <std::vector <Real> > > > w,std::vector <char> aal){
+  this->aa_list=aal;
+  this->naa_=this->aa_list.size();
+  this->aa_dict=init_dict(this->aa_list);
+  this->weights=w;
+}
+
+
+}}}

modules/seq/alg/src/pair_subst_weight_matrix.hh

+//------------------------------------------------------------------------------
+// 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
+//------------------------------------------------------------------------------
+#ifndef OST_SEQ_PAIR_SUBST_WEIGHT_MATRIX_HH
+#define OST_SEQ_PAIR_SUBST_WEIGHT_MATRIX_HH
+
+#include <map>
+#include <ost/seq/alg/module_config.hh>
+
+
+/*
+  Author: Niklaus Johner
+ */
+namespace ost { namespace seq { namespace alg {
+
+  
+/// \brief position-independet pair substitution weight matrix  
+struct DLLEXPORT_OST_SEQ_ALG PairSubstWeightMatrix {
+  std::vector <std::vector <std::vector <std::vector <Real> > > > weights;
+  std::vector <char> aa_list;
+  std::map <char,int> aa_dict;
+  int naa_;
+
+  PairSubstWeightMatrix();
+  PairSubstWeightMatrix(std::vector <std::vector <std::vector <std::vector <Real> > > >,std::vector <char>);
+};
+
+PairSubstWeightMatrix DLLEXPORT_OST_SEQ_ALG LoadDefaultPairSubstWeightMatrix();
+
+}}}
+
+#endif
\ No newline at end of file