Output models in R

From Bio++ Wiki
Jump to: navigation, search


Here is a small code to get the output of generators and equilibrium frequencies for R. For compilation, proceed:

g++ -lbpp-core -lbpp-seq -lbpp-phyl mod2R.cpp -o mod2R

and for usage:

mod2R "alphabet=DNA" "model=GTR(a=3)"

<source lang="cpp">

  1. include <Bpp/App/BppApplication.h>

// The SeqLib library:

  1. include <Bpp/Seq/Alphabet.all>
  2. include <Bpp/Seq/App/SequenceApplicationTools.h>

// The Phyllib library

  1. include <Bpp/Phyl/Model.all>
  2. include <Bpp/Phyl/App/PhylogeneticsApplicationTools.h>


// The STl library:

  1. include <cstdlib>
  2. include <iostream>

using namespace std;

using namespace bpp;

int main(int argc, char* argv[]) {

 VectorSiteContainer* sites=0;
 BppApplication bppml(argc, argv, "BppML");
 Alphabet* alphabet = SequenceApplicationTools::getAlphabet(bppml.getParams(), "", false, false);
 auto_ptr<GeneticCode> gCode;
 CodonAlphabet* codonAlphabet = dynamic_cast<CodonAlphabet*>(alphabet);
 if (codonAlphabet) {
   string codeDesc = ApplicationTools::getStringParameter("genetic_code", bppml.getParams(), "Standard", "", true, true);
   gCode.reset(SequenceApplicationTools::getGeneticCode(codonAlphabet->getNucleicAlphabet(), codeDesc));
 }


 SubstitutionModel* model=PhylogeneticsApplicationTools::getSubstitutionModel(alphabet,gCode.get(),sites,bppml.getParams(), "", true, true);
 
 unsigned int s=alphabet->getSize();
 ParameterList pl=model->getParameters();
 for (unsigned int i=0;i<pl.size();i++)
   cerr << pl[i].getName() << " " << pl[i].getValue() << endl;
 cout << "Alph=c(";
 for (unsigned int j=0;j<s-1;j++)
   cout << "\"" << alphabet->intToChar(j) << "\"" << ",";
 cout << "\"" << alphabet->intToChar(s-1) << "\"" << ")" << endl << endl;
 vector<const SubstitutionModel*> vSM;
 
 if (dynamic_cast<MixedSubstitutionModel*>(model)!=NULL){
   MixedSubstitutionModel* pMSM=dynamic_cast<MixedSubstitutionModel*>(model);
   for (unsigned int i=0;i<pMSM->getNumberOfModels();i++)
     vSM.push_back(pMSM->getNModel(i));
 }
 else
   vSM.push_back(model);
 for (unsigned int nm=0; nm<vSM.size(); nm++)
   {
     const SubstitutionModel* pSM=vSM[nm];
     if (vSM.size()==1)
       cout << "Q=matrix(c(";
     else
       cout << "Q" << nm << "=matrix(c(";
                         
     for (unsigned int i=0;i<s-1;i++){
       for (unsigned int j=0;j<s;j++)
         cout << pSM->Qij(i,j) << ", ";
       cout << endl;
     }
     
     for (unsigned int j=0;j<s-1;j++)
       cout << pSM->Qij(s-1,j) << ", ";
     cout << pSM->Qij(s-1, s-1) << "),ncol=" << s << ",byrow=T)" << endl;
     cout << endl << endl;
     
     const vector<double>& freq=pSM->getFrequencies();
     
     if (vSM.size()==1)
       cout << "freq=c(";
     else
       cout << "freq" << nm << "=c(";
     for (unsigned int j=0;j<s-1;j++)
       cout << freq[j] << ",";
     cout << freq[s-1] << ")" << endl << endl;
   }
 delete alphabet;
 delete model;

}


</source>