Cma.h

Go to the documentation of this file.

/*
 * Cma.h
 *
 * CMA stands for correlated mutation analysis. This class implements
 * information-theoretical measures to compute correlations between
 * positions in a set of aligned sequences.
 * The correlation values computed are Mutual Information (MI) and noiseless
 * Mutual Information (MIp). MIp is a correlation measure proposed in the
 * paper "Mutual information without the influence of phylogeny or entropy
 * dramatically improves residue contact prediction" by S. D. Dunn et al.,
 * Bioinformatics Vol. 24 no. 3, 333-340, 2008.
 * The clustering of correlated positions follows a simple heuristic, which
 * can be found in Thomas Runkler, "Information Mining", chapter 5, 53-58,
 * 2000.
 *
 *  Created on: Dec 16, 2009
 *      Author: Breno Faria
 *
 *  Institute of Microbiology (Technical University Munich)
 *  http://www.arb-home.de/
 */


#ifndef CMA_H
#define CMA_H

#ifndef _GLIBCXX_IOSTREAM
#include <iostream>
#endif
#ifndef _GLIBCXX_VECTOR
#include <vector>
#endif
#ifndef _GLIBCXX_MAP
#include <map>
#endif
#ifndef _GLIBCXX_CMATH
#include <cmath>
#endif
#ifndef _GLIBCXX_STRING
#include <string>
#endif
#ifndef _GLIBCXX_LIST
#include <list>
#endif

#include <eigen/Eigen/Eigen>

#ifndef CXXFORWARD_H
#include <cxxforward.h>
#endif

using namespace std; // @@@ wtf! this really is a nono
using Eigen::VectorXd;
using Eigen::VectorXi;
using Eigen::MatrixXd;

// Vector of maps to store the histogram of sequences.
typedef vector<map<int, int> > MapVecType;
// Vector of vectors to store the aligned sequences (each one represented by a vector)
typedef vector<vector<string> > VecVecType; // @@@ string always is length 1. this should be vector<vector<char>> or vector<string>
// Matrix of maps to store the joint histogram of sequences.
typedef vector<MapVecType> MapMatrixType;

// Struct representing the tuple (pos1, pos2, mutualInformation(pos1,pos2)).
struct MITuple {
    int pos1;
    int pos2;
    double MI;
};

class Cma FINAL_TYPE {

private:
    static CONSTEXPR double epsilon = 1.0e-9;
    static CONSTEXPR double penalty = 5.;
    vector<string> alphabet; // @@@ string always is length 1. this should be vector<char>
    map<string, int> alphabet_map; // @@@ string is length 1 most of the time, disregarding alphabet_map["total"] which is used as global counter. this should be vector<char>
    int num_of_seqs;
    VectorXd entropy;
    MatrixXd JointEntropy;
    MatrixXd MutualInformation;
    MatrixXd MutualInformationP;
    VectorXi clusters;

    void initAlphabet(vector<string> alph);
    MapMatrixType buildJointHistogram(const VecVecType & alignedSequences);
    MatrixXd computeJointEntropy(MapMatrixType & hist);
    MatrixXd computeJointEntropy(const VecVecType & seq);
    VectorXd computeMeanMutualInformationVector();
    double computeOveralMeanMutualInformation();
    bool isValid(string & key1);


public:

    MatrixXd computeMutualInformation(VecVecType & seq);
    MatrixXd computeMutualInformationP(VecVecType & seq);
    VectorXi computeClusters(list<MITuple> mituples, size_t size, double threshold);
    list<MITuple> compute_mituples(MatrixXd mutualInformation);

    MatrixXd getMI();
    MatrixXd getMIp();
    MatrixXd getEntropy();
    MatrixXd getJointEntropy();
    VectorXi getClusters();

    Cma(vector<string> & alph, int seq_number);
    virtual ~Cma();
};

#else
#error Cma.h included twice
#endif // CMA_H