Cma.cxx

Go to the documentation of this file.

/*
 * Cma.cpp
 *
 *  Created on: Dec 16, 2009
 *      Author: Breno Faria
 *
 *  Institute of Microbiology (Technical University Munich)
 *  http://www.arb-home.de/
 */

#include <iostream>
#include <iomanip>
#include <time.h>

#include "Cma.h"

//>--------------------------Global definitions--------------------------------

// Some local function prototypes.
void unifyCluster(int cluster1, int cluster2, VectorXi & result);

//>-----------------------local helper functions-------------------------------

static bool compareMITuples(MITuple first, MITuple second) {
    return first.MI > second.MI;
}

void unifyCluster(int cluster1, int cluster2, VectorXi & result) {
    size_t size = result.size();
    for (size_t i = 0; i < size; ++i) {
        if (result[i] == cluster2) {
            result[i] = cluster1;
        }
    }
}

//>---------------------------Public methods-----------------------------------

Cma::Cma(vector<string> & alph, int seq_num) {

    initAlphabet(alph);
    num_of_seqs = seq_num;

}

Cma::~Cma() {
}

MatrixXd Cma::computeMutualInformation(VecVecType & seq) {

    JointEntropy = computeJointEntropy(seq);

    JointEntropy = JointEntropy + JointEntropy.adjoint() - JointEntropy.diagonal().asDiagonal();

    entropy = JointEntropy.diagonal();

    cout << "computing Mutual Information... ";
    flush(cout);
    clock_t start = clock();

    VectorXd Ones = VectorXd::Ones(entropy.size());
    MatrixXd Hx = (entropy * Ones.transpose()).transpose();
    MatrixXd Hy = Hx.transpose();

    MutualInformation = Hx + Hy - JointEntropy;

    cout << "completed in " << ((double) clock() - start) / CLOCKS_PER_SEC
         << "s." << endl;

    return MutualInformation;

}

MatrixXd Cma::computeMutualInformationP(VecVecType & seq) {

    if (MutualInformation.size() == 0) {
        computeMutualInformation(seq);
    }

    cout << "computing noiseless Mutual Information (MIp)... ";
    flush(cout);
    clock_t start = clock();
    /*
     * We must remove negative MI values (the negative values are artificially
     * generated, see comments in computeJointEntropy).
     */
    for (int i = 0; i < MutualInformation.rows(); ++i) {
        for (int j = 0; j < MutualInformation.cols(); ++j) {
            if (MutualInformation(i, j) < 0.) {
                MutualInformation(i, j) = 0.;
            }
        }
    }

    VectorXd mMIs = computeMeanMutualInformationVector();
    double mMI = computeOveralMeanMutualInformation();

    MutualInformationP = MutualInformation - ((mMIs * mMIs.transpose()) / mMI);

    cout << "completed in " << ((double) clock() - start) / CLOCKS_PER_SEC
         << "s." << endl;

    return MutualInformationP;

}

list<MITuple> Cma::compute_mituples(MatrixXd mutualInformation) {
    size_t size = mutualInformation.cols();

    list<MITuple> mituples;

    // populate the mituples list with all entries in the upper triangular matrix
    // of mutualInformation.
    for (size_t i = 0; i < size; ++i) {
        for (size_t j = i + 1; j < size; ++j) {
            MITuple curr;
            curr.MI = mutualInformation(i, j);
            curr.pos1 = i;
            curr.pos2 = j;
            mituples.push_back(curr);
        }
    }

    // sort mituples by MI-value
    mituples.sort(compareMITuples);

    return mituples;
}

VectorXi Cma::computeClusters(list<MITuple> mituples, size_t size,
                              double threshold) {

    cout << "computing clusters of correlated positions... ";
    flush(cout);

    VectorXi result = VectorXi::Zero(size);

    //compute the clusters.
    int cluster = 1;
    int rest = int(size);
    for (list<MITuple>::iterator it = mituples.begin(); it != mituples.end(); ++it) {
        double result1 = abs(result[it->pos1]);
        double result2 = abs(result[it->pos2]);
        if (it->MI <= threshold || rest == 0) {
            break;
        }
        if (result1 < epsilon and result2 < epsilon) {
            result[it->pos1] = cluster;
            result[it->pos2] = cluster++;
            rest -= 2;
        }
        else if (result1 > epsilon and result2 < epsilon) {
            result[it->pos2] = result[it->pos1];
            rest -= 1;
        }
        else if (result1 < epsilon and result2 > epsilon) {
            result[it->pos1] = result[it->pos2];
            rest -= 1;
        }
        else if (result1 > epsilon and result2 > epsilon && result1 - result2 > epsilon) {
            unifyCluster(result[it->pos1], result[it->pos2], result);
        }
    }

    clusters = result;

    cout << "done." << endl;

    return result;
}

//>-------------------------getters and setters--------------------------------

MatrixXd Cma::getEntropy() {

    return entropy;

}

MatrixXd Cma::getJointEntropy() {

    return JointEntropy;

}

MatrixXd Cma::getMI() {

    return MutualInformation;

}

MatrixXd Cma::getMIp() {

    return MutualInformationP;

}

VectorXi Cma::getClusters() {

    return clusters;

}

//>----------------------------private methods---------------------------------

void Cma::initAlphabet(vector<string> alph) {

    alphabet = vector<string> (alph);

    int i = 0;
    for (vector<string>::iterator it = alph.begin(); it != alph.end(); ++it) {
        alphabet_map[*it] = i;
        i++;
    }
    for (vector<string>::iterator it = alph.begin(); it != alph.end(); ++it) {
        for (vector<string>::iterator it2 = alph.begin(); it2 != alph.end(); ++it2) {
            alphabet_map[*it + *it2] = i;
            i++;
        }
    }
    alphabet_map["total"] = i;
}

MatrixXd Cma::computeJointEntropy(const VecVecType & seqs) {

    MapMatrixType hist = buildJointHistogram(seqs);
    return Cma::computeJointEntropy(hist);

}

MatrixXd Cma::computeJointEntropy(MapMatrixType & hist) {

    cout << "computing joint entropy (this may take a while)... ";
    flush(cout);
    clock_t start = clock();

    int hist_size = int(hist.size());
    MatrixXd result(hist_size, hist_size);
    result.setZero(hist_size, hist_size);
    for (unsigned int i = 0; i < (unsigned int) hist_size; ++i) {

        //progress "bar"
        if (hist_size > 30 and i % (hist_size / 30) == 0) {
            cout << "(" << setw(6) << setiosflags(ios::fixed)
                 << setprecision(2) << i / float(hist_size) * 100 << "%)";
            flush(cout);
        }

        for (unsigned int j = i; j < hist[i].size(); ++j) {
            int total_i_j = hist[i][j][alphabet_map.at("total")];
            int total_i_i = hist[i][i][alphabet_map.at("total")];
            int total_j_j = hist[j][j][alphabet_map.at("total")];
            int mismatches = total_i_i + total_j_j - 2 * total_i_j;
            if (total_i_j != 0) {
                double result_i_j = 0.;
                for (map<int, int>::iterator it = hist[i][j].begin(); it
                         != hist[i][j].end(); ++it) {
                    double pair = double(it->second) + Cma::epsilon;
                    result_i_j += (pair / total_i_j) * log2(pair / total_i_j);
                }

                result(i, j) = -result_i_j + double(mismatches) / (total_i_i
                                                                   + total_j_j) * Cma::penalty;

            } else {
                result(i, j) = 5.;
            }

        }
        //progress "bar"
        if (hist_size > 30 and i % (hist_size / 30) == 0) {
            cout << "\b\b\b\b\b\b\b\b\b";
        }
    }
    flush(cout);
    cout << "completed in " << ((double) clock() - start) / CLOCKS_PER_SEC
         << "s." << endl;
    return result;

}

bool Cma::isValid(string & key) {

    map<string, int>::iterator it = alphabet_map.find(key);
    return it != alphabet_map.end();

}

MapMatrixType Cma::buildJointHistogram(const VecVecType & alignedSequences) {

    cout << "building histogram (this may take a while)... ";
    flush(cout);
    clock_t start = clock();

    if (alignedSequences.empty()) {
        throw "no (marked?) sequences";
    }

    // here we defined one dimension of the matrix as being the same as the
    // length of the sequences.
    MapMatrixType result(alignedSequences[0].size(), MapVecType(
                             alignedSequences[0].size()));

    int ii = 0;
    // for each sequence:
    for (VecVecType::const_iterator it = alignedSequences.begin(); it
             != alignedSequences.end(); ++it) {

        cout << "(" << setw(6) << setiosflags(ios::fixed) << setprecision(2)
             << ii / float(num_of_seqs) * 100 << "%)";
        flush(cout);

        // for each first position in the sequence
        for (size_t i = 0; i < it->size(); ++i) {
            string key1 = it->at(i);
            if (isValid(key1)) {

                //for each second position in the sequence
                for (size_t j = i; j < it->size(); ++j) {
                    string key2 = it->at(j);
                    if (isValid(key2)) {
                        result[i][j][alphabet_map.at(key1 + key2)] += 1;
                        result[i][j][alphabet_map.at("total")] += 1;
                    }
                }
            }
        }
        cout << "\b\b\b\b\b\b\b\b\b";
        ii++;
    }
    flush(cout);
    cout << "completed in " << ((double) clock() - start) / CLOCKS_PER_SEC
         << "s." << endl;

    return result;
}

VectorXd Cma::computeMeanMutualInformationVector() {

    return MutualInformation.rowwise().sum() / num_of_seqs;

}

double Cma::computeOveralMeanMutualInformation() {

    return 2. * MutualInformation.sum() / num_of_seqs / num_of_seqs;

}