binaries/src/mafft/extensions/mxscarna_src/BPPMatrix.hpp

   1 //////////////////////////////////////////////////////////////
   2 // BPPMatrix.hpp
   3 //
   4 // save the Base Pairing Probability Matrix for each sequences
   5 //////////////////////////////////////////////////////////////
   6
   7 #ifndef __BBPMatrix_HPP__
   8 #define __BBPMatrix_HPP__
   9
  10 #include <iostream>
  11 #include <string>
  12 #include "SparseMatrix.h"
  13 #include "McCaskill.hpp"
  14 #include "nrutil.h"
  15
  16
  17
  18 using namespace std;
  19
  20 namespace MXSCARNA{
  21 class BPPMatrix {
  22 private:
  23
  24     int seqLength;       // sequence length;
  25     float cutOff;        // the threshold of probability
  26
  27     BPPMatrix() {}
  28 public:
  29     SparseMatrix bppMat; // base pairing probability matrix 1-origin(1..seqLength)
  30     BPPMatrix(int bppmode, const string &seq, int seqLength, float inCutOff = 0.0001) : seqLength (seqLength), cutOff(inCutOff) {
  31       if( bppmode == 'r' ) // by katoh
  32         readBPPMatrix(seq);
  33       else if( bppmode == 'w' )
  34         setandwriteBPPMatrix(seq);
  35       else
  36         setBPPMatrix(seq);
  37     }
  38     BPPMatrix(int seqLength, float inCutOff, const Trimat<float> & tmpBppMat) : seqLength(seqLength), cutOff(inCutOff) {
  39       bppMat.SetSparseMatrix(seqLength, seqLength, tmpBppMat, cutOff);
  40     }
  41
  42
  43     void setBPPMatrix(const string &seq) {
  44         const char *tmpSeq = seq.c_str();
  45         McCaskill mc(seqLength, &tmpSeq[1]);
  46         mc.calcPartitionFunction();
  47         Trimat<float> tmpBppMat(seqLength + 1);
  48
  49         for (int i = 0; i < seqLength; i++) {
  50           for(int j = i; j < seqLength; j++) {
  51             tmpBppMat.ref(i+1, j+1) = mc.getProb(i,j);
  52           }
  53         }
  54         bppMat.SetSparseMatrix(seqLength, seqLength, tmpBppMat, cutOff);
  55     }
  56
  57     void setandwriteBPPMatrix(const string &seq) { // by katoh
  58         const char *tmpSeq = seq.c_str();
  59         McCaskill mc(seqLength, &tmpSeq[1]);
  60         mc.calcPartitionFunction();
  61         Trimat<float> tmpBppMat(seqLength + 1);
  62         float tmpbp;
  63
  64         fprintf( stdout, ">\n" );
  65         for (int i = 0; i < seqLength; i++) {
  66           for(int j = i; j < seqLength; j++) {
  67         tmpbp = mc.getProb(i,j);
  68         if (tmpbp > 0.0001)
  69           fprintf( stdout, "%d %d %50.40f\n", i, j, tmpbp );
  70             tmpBppMat.ref(i+1, j+1) = tmpbp;
  71           }
  72         }
  73         bppMat.SetSparseMatrix(seqLength, seqLength, tmpBppMat, cutOff);
  74     }
  75
  76     void readBPPMatrix(const string &seq) { // by katoh
  77         const char *tmpSeq = seq.c_str();
  78         char oneline[1000];
  79         int onechar;
  80         float prob;
  81         int posi, posj;
  82         static FILE *bppfp = NULL;
  83
  84
  85         if( bppfp == NULL )
  86         {
  87                 bppfp = fopen( "_bpp", "r" );
  88                 if( bppfp == NULL )
  89                 {
  90                         fprintf( stderr, "Cannot open _bpp.\n" );
  91                         exit( 1 );
  92                 }
  93         }
  94
  95         Trimat<float> tmpBppMat(seqLength + 1);
  96         fgets( oneline, 999, bppfp );
  97         if( oneline[0] != '>' )
  98         {
  99                 fprintf( stderr, "Format error\n" );
 100                 exit( 1 );
 101         }
 102         while( 1 )
 103         {
 104                 onechar = getc( bppfp );
 105                 ungetc( onechar, bppfp );
 106                 if( onechar == '>' || onechar == EOF ) break;
 107
 108                 fgets( oneline, 999, bppfp );
 109                 sscanf( oneline, "%d %d %f", &posi, &posj, &prob );
 110 //              fprintf( stderr, "%d %d -> %f\n", posi, posj, prob );
 111                 tmpBppMat.ref(posi+1, posj+1) = prob;
 112         }
 113
 114         bppMat.SetSparseMatrix(seqLength, seqLength, tmpBppMat, cutOff);
 115     }
 116
 117     float GetProb(int i, int j) {
 118         return bppMat.GetValue(i,j);
 119     }
 120
 121     float GetLength() const {
 122         return seqLength;
 123     }
 124
 125     void updateBPPMatrix(const Trimat<float> &inbppMat) {
 126         bppMat.SetSparseMatrix(seqLength, seqLength, inbppMat, cutOff);
 127     }
 128 };
 129 }
 130 #endif // __BPPMatrix_HPP__