--- /dev/null
+/////////////////////////////////////////////////////////////////
+// ProbabilisticModel.h
+//
+// Routines for (1) posterior probability computations
+// (2) chained anchoring
+// (3) maximum weight trace alignment
+/////////////////////////////////////////////////////////////////
+
+#ifndef PROBABILISTICMODEL_H
+#define PROBABILISTICMODEL_H
+
+#include <list>
+#include <cmath>
+#include <cstdio>
+#include "SafeVector.h"
+#include "ScoreType.h"
+#include "SparseMatrix.h"
+#include "MultiSequence.h"
+
+using namespace std;
+
+const int NumMatchStates = 1; // note that in this version the number
+ // of match states is fixed at 1...will
+ // change in future versions
+const int NumMatrixTypes = NumMatchStates + NumInsertStates * 2;
+
+/////////////////////////////////////////////////////////////////
+// ProbabilisticModel
+//
+// Class for storing the parameters of a probabilistic model and
+// performing different computations based on those parameters.
+// In particular, this class handles the computation of
+// posterior probabilities that may be used in alignment.
+/////////////////////////////////////////////////////////////////
+
+class ProbabilisticModel {
+
+ float initialDistribution[NumMatrixTypes]; // holds the initial probabilities for each state
+ float transProb[NumMatrixTypes][NumMatrixTypes]; // holds all state-to-state transition probabilities
+ float matchProb[256][256]; // emission probabilities for match states
+ float insProb[256][NumMatrixTypes]; // emission probabilities for insert states
+
+ public:
+
+ /////////////////////////////////////////////////////////////////
+ // ProbabilisticModel::ProbabilisticModel()
+ //
+ // Constructor. Builds a new probabilistic model using the
+ // given parameters.
+ /////////////////////////////////////////////////////////////////
+
+ ProbabilisticModel (const VF &initDistribMat, const VF &gapOpen, const VF &gapExtend,
+ const VVF &emitPairs, const VF &emitSingle){
+
+ // build transition matrix
+ VVF transMat (NumMatrixTypes, VF (NumMatrixTypes, 0.0f));
+ transMat[0][0] = 1;
+ for (int i = 0; i < NumInsertStates; i++){
+ transMat[0][2*i+1] = gapOpen[2*i];
+ transMat[0][2*i+2] = gapOpen[2*i+1];
+ transMat[0][0] -= (gapOpen[2*i] + gapOpen[2*i+1]);
+ assert (transMat[0][0] > 0);
+ transMat[2*i+1][2*i+1] = gapExtend[2*i];
+ transMat[2*i+2][2*i+2] = gapExtend[2*i+1];
+ transMat[2*i+1][2*i+2] = 0;
+ transMat[2*i+2][2*i+1] = 0;
+ transMat[2*i+1][0] = 1 - gapExtend[2*i];
+ transMat[2*i+2][0] = 1 - gapExtend[2*i+1];
+ }
+
+ // create initial and transition probability matrices
+ for (int i = 0; i < NumMatrixTypes; i++){
+ initialDistribution[i] = LOG (initDistribMat[i]);
+ for (int j = 0; j < NumMatrixTypes; j++)
+ transProb[i][j] = LOG (transMat[i][j]);
+ }
+
+ // create insertion and match probability matrices
+ for (int i = 0; i < 256; i++){
+ for (int j = 0; j < NumMatrixTypes; j++)
+ insProb[i][j] = LOG (emitSingle[i]);
+ for (int j = 0; j < 256; j++)
+ matchProb[i][j] = LOG (emitPairs[i][j]);
+ }
+ }
+
+ /////////////////////////////////////////////////////////////////
+ // ProbabilisticModel::ComputeForwardMatrix()
+ //
+ // Computes a set of forward probability matrices for aligning
+ // seq1 and seq2.
+ //
+ // For efficiency reasons, a single-dimensional floating-point
+ // array is used here, with the following indexing scheme:
+ //
+ // forward[i + NumMatrixTypes * (j * (seq2Length+1) + k)]
+ // refers to the probability of aligning through j characters
+ // of the first sequence, k characters of the second sequence,
+ // and ending in state i.
+ /////////////////////////////////////////////////////////////////
+
+ VF *ComputeForwardMatrix (Sequence *seq1, Sequence *seq2) const {
+
+ assert (seq1);
+ assert (seq2);
+
+ const int seq1Length = seq1->GetLength();
+ const int seq2Length = seq2->GetLength();
+
+ // retrieve the points to the beginning of each sequence
+ SafeVector<char>::iterator iter1 = seq1->GetDataPtr();
+ SafeVector<char>::iterator iter2 = seq2->GetDataPtr();
+
+ // create matrix
+ VF *forwardPtr = new VF (NumMatrixTypes * (seq1Length+1) * (seq2Length+1), LOG_ZERO);
+ assert (forwardPtr);
+ VF &forward = *forwardPtr;
+
+ // initialization condition
+ forward[0 + NumMatrixTypes * (1 * (seq2Length+1) + 1)] =
+ initialDistribution[0] + matchProb[(unsigned char) iter1[1]][(unsigned char) iter2[1]];
+
+ for (int k = 0; k < NumInsertStates; k++){
+ forward[2*k+1 + NumMatrixTypes * (1 * (seq2Length+1) + 0)] =
+ initialDistribution[2*k+1] + insProb[(unsigned char) iter1[1]][k];
+ forward[2*k+2 + NumMatrixTypes * (0 * (seq2Length+1) + 1)] =
+ initialDistribution[2*k+2] + insProb[(unsigned char) iter2[1]][k];
+ }
+
+ // remember offset for each index combination
+ int ij = 0;
+ int i1j = -seq2Length - 1;
+ int ij1 = -1;
+ int i1j1 = -seq2Length - 2;
+
+ ij *= NumMatrixTypes;
+ i1j *= NumMatrixTypes;
+ ij1 *= NumMatrixTypes;
+ i1j1 *= NumMatrixTypes;
+
+ // compute forward scores
+ for (int i = 0; i <= seq1Length; i++){
+ unsigned char c1 = (i == 0) ? '~' : (unsigned char) iter1[i];
+ for (int j = 0; j <= seq2Length; j++){
+ unsigned char c2 = (j == 0) ? '~' : (unsigned char) iter2[j];
+
+ if (i > 1 || j > 1){
+ if (i > 0 && j > 0){
+ forward[0 + ij] = forward[0 + i1j1] + transProb[0][0];
+ for (int k = 1; k < NumMatrixTypes; k++)
+ LOG_PLUS_EQUALS (forward[0 + ij], forward[k + i1j1] + transProb[k][0]);
+ forward[0 + ij] += matchProb[c1][c2];
+ }
+ if (i > 0){
+ for (int k = 0; k < NumInsertStates; k++)
+ forward[2*k+1 + ij] = insProb[c1][k] +
+ LOG_ADD (forward[0 + i1j] + transProb[0][2*k+1],
+ forward[2*k+1 + i1j] + transProb[2*k+1][2*k+1]);
+ }
+ if (j > 0){
+ for (int k = 0; k < NumInsertStates; k++)
+ forward[2*k+2 + ij] = insProb[c2][k] +
+ LOG_ADD (forward[0 + ij1] + transProb[0][2*k+2],
+ forward[2*k+2 + ij1] + transProb[2*k+2][2*k+2]);
+ }
+ }
+
+ ij += NumMatrixTypes;
+ i1j += NumMatrixTypes;
+ ij1 += NumMatrixTypes;
+ i1j1 += NumMatrixTypes;
+ }
+ }
+
+ return forwardPtr;
+ }
+
+ /////////////////////////////////////////////////////////////////
+ // ProbabilisticModel::ComputeBackwardMatrix()
+ //
+ // Computes a set of backward probability matrices for aligning
+ // seq1 and seq2.
+ //
+ // For efficiency reasons, a single-dimensional floating-point
+ // array is used here, with the following indexing scheme:
+ //
+ // backward[i + NumMatrixTypes * (j * (seq2Length+1) + k)]
+ // refers to the probability of starting in state i and
+ // aligning from character j+1 to the end of the first
+ // sequence and from character k+1 to the end of the second
+ // sequence.
+ /////////////////////////////////////////////////////////////////
+
+ VF *ComputeBackwardMatrix (Sequence *seq1, Sequence *seq2) const {
+
+ assert (seq1);
+ assert (seq2);
+
+ const int seq1Length = seq1->GetLength();
+ const int seq2Length = seq2->GetLength();
+ SafeVector<char>::iterator iter1 = seq1->GetDataPtr();
+ SafeVector<char>::iterator iter2 = seq2->GetDataPtr();
+
+ // create matrix
+ VF *backwardPtr = new VF (NumMatrixTypes * (seq1Length+1) * (seq2Length+1), LOG_ZERO);
+ assert (backwardPtr);
+ VF &backward = *backwardPtr;
+
+ // initialization condition
+ for (int k = 0; k < NumMatrixTypes; k++)
+ backward[NumMatrixTypes * ((seq1Length+1) * (seq2Length+1) - 1) + k] = initialDistribution[k];
+
+ // remember offset for each index combination
+ int ij = (seq1Length+1) * (seq2Length+1) - 1;
+ int i1j = ij + seq2Length + 1;
+ int ij1 = ij + 1;
+ int i1j1 = ij + seq2Length + 2;
+
+ ij *= NumMatrixTypes;
+ i1j *= NumMatrixTypes;
+ ij1 *= NumMatrixTypes;
+ i1j1 *= NumMatrixTypes;
+
+ // compute backward scores
+ for (int i = seq1Length; i >= 0; i--){
+ unsigned char c1 = (i == seq1Length) ? '~' : (unsigned char) iter1[i+1];
+ for (int j = seq2Length; j >= 0; j--){
+ unsigned char c2 = (j == seq2Length) ? '~' : (unsigned char) iter2[j+1];
+
+ if (i < seq1Length && j < seq2Length){
+ const float ProbXY = backward[0 + i1j1] + matchProb[c1][c2];
+ for (int k = 0; k < NumMatrixTypes; k++)
+ LOG_PLUS_EQUALS (backward[k + ij], ProbXY + transProb[k][0]);
+ }
+ if (i < seq1Length){
+ for (int k = 0; k < NumInsertStates; k++){
+ LOG_PLUS_EQUALS (backward[0 + ij], backward[2*k+1 + i1j] + insProb[c1][k] + transProb[0][2*k+1]);
+ LOG_PLUS_EQUALS (backward[2*k+1 + ij], backward[2*k+1 + i1j] + insProb[c1][k] + transProb[2*k+1][2*k+1]);
+ }
+ }
+ if (j < seq2Length){
+ for (int k = 0; k < NumInsertStates; k++){
+ LOG_PLUS_EQUALS (backward[0 + ij], backward[2*k+2 + ij1] + insProb[c2][k] + transProb[0][2*k+2]);
+ LOG_PLUS_EQUALS (backward[2*k+2 + ij], backward[2*k+2 + ij1] + insProb[c2][k] + transProb[2*k+2][2*k+2]);
+ }
+ }
+
+ ij -= NumMatrixTypes;
+ i1j -= NumMatrixTypes;
+ ij1 -= NumMatrixTypes;
+ i1j1 -= NumMatrixTypes;
+ }
+ }
+
+ return backwardPtr;
+ }
+
+ /////////////////////////////////////////////////////////////////
+ // ProbabilisticModel::ComputeTotalProbability()
+ //
+ // Computes the total probability of an alignment given
+ // the forward and backward matrices.
+ /////////////////////////////////////////////////////////////////
+
+ float ComputeTotalProbability (int seq1Length, int seq2Length,
+ const VF &forward, const VF &backward) const {
+
+ // compute total probability
+ float totalForwardProb = LOG_ZERO;
+ float totalBackwardProb = LOG_ZERO;
+ for (int k = 0; k < NumMatrixTypes; k++){
+ LOG_PLUS_EQUALS (totalForwardProb,
+ forward[k + NumMatrixTypes * ((seq1Length+1) * (seq2Length+1) - 1)] +
+ backward[k + NumMatrixTypes * ((seq1Length+1) * (seq2Length+1) - 1)]);
+ }
+
+ totalBackwardProb =
+ forward[0 + NumMatrixTypes * (1 * (seq2Length+1) + 1)] +
+ backward[0 + NumMatrixTypes * (1 * (seq2Length+1) + 1)];
+
+ for (int k = 0; k < NumInsertStates; k++){
+ LOG_PLUS_EQUALS (totalBackwardProb,
+ forward[2*k+1 + NumMatrixTypes * (1 * (seq2Length+1) + 0)] +
+ backward[2*k+1 + NumMatrixTypes * (1 * (seq2Length+1) + 0)]);
+ LOG_PLUS_EQUALS (totalBackwardProb,
+ forward[2*k+2 + NumMatrixTypes * (0 * (seq2Length+1) + 1)] +
+ backward[2*k+2 + NumMatrixTypes * (0 * (seq2Length+1) + 1)]);
+ }
+
+ // cerr << totalForwardProb << " " << totalBackwardProb << endl;
+
+ return (totalForwardProb + totalBackwardProb) / 2;
+ }
+
+ /////////////////////////////////////////////////////////////////
+ // ProbabilisticModel::ComputePosteriorMatrix()
+ //
+ // Computes the posterior probability matrix based on
+ // the forward and backward matrices.
+ /////////////////////////////////////////////////////////////////
+
+ VF *ComputePosteriorMatrix (Sequence *seq1, Sequence *seq2,
+ const VF &forward, const VF &backward) const {
+
+ assert (seq1);
+ assert (seq2);
+
+ const int seq1Length = seq1->GetLength();
+ const int seq2Length = seq2->GetLength();
+
+ float totalProb = ComputeTotalProbability (seq1Length, seq2Length,
+ forward, backward);
+
+ // compute posterior matrices
+ VF *posteriorPtr = new VF((seq1Length+1) * (seq2Length+1)); assert (posteriorPtr);
+ VF &posterior = *posteriorPtr;
+
+ int ij = 0;
+ VF::iterator ptr = posterior.begin();
+
+ for (int i = 0; i <= seq1Length; i++){
+ for (int j = 0; j <= seq2Length; j++){
+ *(ptr++) = EXP (min (LOG_ONE, forward[ij] + backward[ij] - totalProb));
+ ij += NumMatrixTypes;
+ }
+ }
+
+ posterior[0] = 0;
+
+ return posteriorPtr;
+ }
+
+ /*
+ /////////////////////////////////////////////////////////////////
+ // ProbabilisticModel::ComputeExpectedCounts()
+ //
+ // Computes the expected counts for the various transitions.
+ /////////////////////////////////////////////////////////////////
+
+ VVF *ComputeExpectedCounts () const {
+
+ assert (seq1);
+ assert (seq2);
+
+ const int seq1Length = seq1->GetLength();
+ const int seq2Length = seq2->GetLength();
+ SafeVector<char>::iterator iter1 = seq1->GetDataPtr();
+ SafeVector<char>::iterator iter2 = seq2->GetDataPtr();
+
+ // compute total probability
+ float totalProb = ComputeTotalProbability (seq1Length, seq2Length,
+ forward, backward);
+
+ // initialize expected counts
+ VVF *countsPtr = new VVF(NumMatrixTypes + 1, VF(NumMatrixTypes, LOG_ZERO)); assert (countsPtr);
+ VVF &counts = *countsPtr;
+
+ // remember offset for each index combination
+ int ij = 0;
+ int i1j = -seq2Length - 1;
+ int ij1 = -1;
+ int i1j1 = -seq2Length - 2;
+
+ ij *= NumMatrixTypes;
+ i1j *= NumMatrixTypes;
+ ij1 *= NumMatrixTypes;
+ i1j1 *= NumMatrixTypes;
+
+ // compute expected counts
+ for (int i = 0; i <= seq1Length; i++){
+ unsigned char c1 = (i == 0) ? '~' : (unsigned char) iter1[i];
+ for (int j = 0; j <= seq2Length; j++){
+ unsigned char c2 = (j == 0) ? '~' : (unsigned char) iter2[j];
+
+ if (i > 0 && j > 0){
+ for (int k = 0; k < NumMatrixTypes; k++)
+ LOG_PLUS_EQUALS (counts[k][0],
+ forward[k + i1j1] + transProb[k][0] +
+ matchProb[c1][c2] + backward[0 + ij]);
+ }
+ if (i > 0){
+ for (int k = 0; k < NumInsertStates; k++){
+ LOG_PLUS_EQUALS (counts[0][2*k+1],
+ forward[0 + i1j] + transProb[0][2*k+1] +
+ insProb[c1][k] + backward[2*k+1 + ij]);
+ LOG_PLUS_EQUALS (counts[2*k+1][2*k+1],
+ forward[2*k+1 + i1j] + transProb[2*k+1][2*k+1] +
+ insProb[c1][k] + backward[2*k+1 + ij]);
+ }
+ }
+ if (j > 0){
+ for (int k = 0; k < NumInsertStates; k++){
+ LOG_PLUS_EQUALS (counts[0][2*k+2],
+ forward[0 + ij1] + transProb[0][2*k+2] +
+ insProb[c2][k] + backward[2*k+2 + ij]);
+ LOG_PLUS_EQUALS (counts[2*k+2][2*k+2],
+ forward[2*k+2 + ij1] + transProb[2*k+2][2*k+2] +
+ insProb[c2][k] + backward[2*k+2 + ij]);
+ }
+ }
+
+ ij += NumMatrixTypes;
+ i1j += NumMatrixTypes;
+ ij1 += NumMatrixTypes;
+ i1j1 += NumMatrixTypes;
+ }
+ }
+
+ // scale all expected counts appropriately
+ for (int i = 0; i < NumMatrixTypes; i++)
+ for (int j = 0; j < NumMatrixTypes; j++)
+ counts[i][j] -= totalProb;
+
+ }
+ */
+
+ /////////////////////////////////////////////////////////////////
+ // ProbabilisticModel::ComputeNewParameters()
+ //
+ // Computes a new parameter set based on the expected counts
+ // given.
+ /////////////////////////////////////////////////////////////////
+
+ void ComputeNewParameters (Sequence *seq1, Sequence *seq2,
+ const VF &forward, const VF &backward,
+ VF &initDistribMat, VF &gapOpen,
+ VF &gapExtend, VVF &emitPairs, VF &emitSingle, bool enableTrainEmissions) const {
+
+ assert (seq1);
+ assert (seq2);
+
+ const int seq1Length = seq1->GetLength();
+ const int seq2Length = seq2->GetLength();
+ SafeVector<char>::iterator iter1 = seq1->GetDataPtr();
+ SafeVector<char>::iterator iter2 = seq2->GetDataPtr();
+
+ // compute total probability
+ float totalProb = ComputeTotalProbability (seq1Length, seq2Length,
+ forward, backward);
+
+ // initialize expected counts
+ VVF transCounts (NumMatrixTypes, VF (NumMatrixTypes, LOG_ZERO));
+ VF initCounts (NumMatrixTypes, LOG_ZERO);
+ VVF pairCounts (256, VF (256, LOG_ZERO));
+ VF singleCounts (256, LOG_ZERO);
+
+ // remember offset for each index combination
+ int ij = 0;
+ int i1j = -seq2Length - 1;
+ int ij1 = -1;
+ int i1j1 = -seq2Length - 2;
+
+ ij *= NumMatrixTypes;
+ i1j *= NumMatrixTypes;
+ ij1 *= NumMatrixTypes;
+ i1j1 *= NumMatrixTypes;
+
+ // compute initial distribution posteriors
+ initCounts[0] = LOG_ADD (forward[0 + NumMatrixTypes * (1 * (seq2Length+1) + 1)] +
+ backward[0 + NumMatrixTypes * (1 * (seq2Length+1) + 1)],
+ forward[0 + NumMatrixTypes * ((seq1Length+1) * (seq2Length+1) - 1)] +
+ backward[0 + NumMatrixTypes * ((seq1Length+1) * (seq2Length+1) - 1)]);
+ for (int k = 0; k < NumInsertStates; k++){
+ initCounts[2*k+1] = LOG_ADD (forward[2*k+1 + NumMatrixTypes * (1 * (seq2Length+1) + 0)] +
+ backward[2*k+1 + NumMatrixTypes * (1 * (seq2Length+1) + 0)],
+ forward[2*k+1 + NumMatrixTypes * ((seq1Length+1) * (seq2Length+1) - 1)] +
+ backward[2*k+1 + NumMatrixTypes * ((seq1Length+1) * (seq2Length+1) - 1)]);
+ initCounts[2*k+2] = LOG_ADD (forward[2*k+2 + NumMatrixTypes * (0 * (seq2Length+1) + 1)] +
+ backward[2*k+2 + NumMatrixTypes * (0 * (seq2Length+1) + 1)],
+ forward[2*k+2 + NumMatrixTypes * ((seq1Length+1) * (seq2Length+1) - 1)] +
+ backward[2*k+2 + NumMatrixTypes * ((seq1Length+1) * (seq2Length+1) - 1)]);
+ }
+
+ // compute expected counts
+ for (int i = 0; i <= seq1Length; i++){
+ unsigned char c1 = (i == 0) ? '~' : (unsigned char) toupper(iter1[i]);
+ for (int j = 0; j <= seq2Length; j++){
+ unsigned char c2 = (j == 0) ? '~' : (unsigned char) toupper(iter2[j]);
+
+ if (i > 0 && j > 0){
+ if (enableTrainEmissions && i == 1 && j == 1){
+ LOG_PLUS_EQUALS (pairCounts[c1][c2],
+ initialDistribution[0] + matchProb[c1][c2] + backward[0 + ij]);
+ LOG_PLUS_EQUALS (pairCounts[c2][c1],
+ initialDistribution[0] + matchProb[c2][c1] + backward[0 + ij]);
+ }
+
+ for (int k = 0; k < NumMatrixTypes; k++){
+ LOG_PLUS_EQUALS (transCounts[k][0],
+ forward[k + i1j1] + transProb[k][0] +
+ matchProb[c1][c2] + backward[0 + ij]);
+ if (enableTrainEmissions && i != 1 || j != 1){
+ LOG_PLUS_EQUALS (pairCounts[c1][c2],
+ forward[k + i1j1] + transProb[k][0] +
+ matchProb[c1][c2] + backward[0 + ij]);
+ LOG_PLUS_EQUALS (pairCounts[c2][c1],
+ forward[k + i1j1] + transProb[k][0] +
+ matchProb[c2][c1] + backward[0 + ij]);
+ }
+ }
+ }
+ if (i > 0){
+ for (int k = 0; k < NumInsertStates; k++){
+ LOG_PLUS_EQUALS (transCounts[0][2*k+1],
+ forward[0 + i1j] + transProb[0][2*k+1] +
+ insProb[c1][k] + backward[2*k+1 + ij]);
+ LOG_PLUS_EQUALS (transCounts[2*k+1][2*k+1],
+ forward[2*k+1 + i1j] + transProb[2*k+1][2*k+1] +
+ insProb[c1][k] + backward[2*k+1 + ij]);
+ if (enableTrainEmissions){
+ if (i == 1 && j == 0){
+ LOG_PLUS_EQUALS (singleCounts[c1],
+ initialDistribution[2*k+1] + insProb[c1][k] + backward[2*k+1 + ij]);
+ }
+ else {
+ LOG_PLUS_EQUALS (singleCounts[c1],
+ forward[0 + i1j] + transProb[0][2*k+1] +
+ insProb[c1][k] + backward[2*k+1 + ij]);
+ LOG_PLUS_EQUALS (singleCounts[c1],
+ forward[2*k+1 + i1j] + transProb[2*k+1][2*k+1] +
+ insProb[c1][k] + backward[2*k+1 + ij]);
+ }
+ }
+ }
+ }
+ if (j > 0){
+ for (int k = 0; k < NumInsertStates; k++){
+ LOG_PLUS_EQUALS (transCounts[0][2*k+2],
+ forward[0 + ij1] + transProb[0][2*k+2] +
+ insProb[c2][k] + backward[2*k+2 + ij]);
+ LOG_PLUS_EQUALS (transCounts[2*k+2][2*k+2],
+ forward[2*k+2 + ij1] + transProb[2*k+2][2*k+2] +
+ insProb[c2][k] + backward[2*k+2 + ij]);
+ if (enableTrainEmissions){
+ if (i == 0 && j == 1){
+ LOG_PLUS_EQUALS (singleCounts[c2],
+ initialDistribution[2*k+2] + insProb[c2][k] + backward[2*k+2 + ij]);
+ }
+ else {
+ LOG_PLUS_EQUALS (singleCounts[c2],
+ forward[0 + ij1] + transProb[0][2*k+2] +
+ insProb[c2][k] + backward[2*k+2 + ij]);
+ LOG_PLUS_EQUALS (singleCounts[c2],
+ forward[2*k+2 + ij1] + transProb[2*k+2][2*k+2] +
+ insProb[c2][k] + backward[2*k+2 + ij]);
+ }
+ }
+ }
+ }
+
+ ij += NumMatrixTypes;
+ i1j += NumMatrixTypes;
+ ij1 += NumMatrixTypes;
+ i1j1 += NumMatrixTypes;
+ }
+ }
+
+ // scale all expected counts appropriately
+ for (int i = 0; i < NumMatrixTypes; i++){
+ initCounts[i] -= totalProb;
+ for (int j = 0; j < NumMatrixTypes; j++)
+ transCounts[i][j] -= totalProb;
+ }
+ if (enableTrainEmissions){
+ for (int i = 0; i < 256; i++){
+ for (int j = 0; j < 256; j++)
+ pairCounts[i][j] -= totalProb;
+ singleCounts[i] -= totalProb;
+ }
+ }
+
+ // compute new initial distribution
+ float totalInitDistribCounts = 0;
+ for (int i = 0; i < NumMatrixTypes; i++)
+ totalInitDistribCounts += exp (initCounts[i]); // should be 2
+ initDistribMat[0] = min (1.0f, max (0.0f, (float) exp (initCounts[0]) / totalInitDistribCounts));
+ for (int k = 0; k < NumInsertStates; k++){
+ float val = (exp (initCounts[2*k+1]) + exp (initCounts[2*k+2])) / 2;
+ initDistribMat[2*k+1] = initDistribMat[2*k+2] = min (1.0f, max (0.0f, val / totalInitDistribCounts));
+ }
+
+ // compute total counts for match state
+ float inMatchStateCounts = 0;
+ for (int i = 0; i < NumMatrixTypes; i++)
+ inMatchStateCounts += exp (transCounts[0][i]);
+ for (int i = 0; i < NumInsertStates; i++){
+
+ // compute total counts for gap state
+ float inGapStateCounts =
+ exp (transCounts[2*i+1][0]) +
+ exp (transCounts[2*i+1][2*i+1]) +
+ exp (transCounts[2*i+2][0]) +
+ exp (transCounts[2*i+2][2*i+2]);
+
+ gapOpen[2*i] = gapOpen[2*i+1] =
+ (exp (transCounts[0][2*i+1]) +
+ exp (transCounts[0][2*i+2])) /
+ (2 * inMatchStateCounts);
+
+ gapExtend[2*i] = gapExtend[2*i+1] =
+ (exp (transCounts[2*i+1][2*i+1]) +
+ exp (transCounts[2*i+2][2*i+2])) /
+ inGapStateCounts;
+ }
+
+ if (enableTrainEmissions){
+ float totalPairCounts = 0;
+ float totalSingleCounts = 0;
+ for (int i = 0; i < 256; i++){
+ for (int j = 0; j <= i; j++)
+ totalPairCounts += exp (pairCounts[j][i]);
+ totalSingleCounts += exp (singleCounts[i]);
+ }
+
+ for (int i = 0; i < 256; i++) if (!islower ((char) i)){
+ int li = (int)((unsigned char) tolower ((char) i));
+ for (int j = 0; j <= i; j++) if (!islower ((char) j)){
+ int lj = (int)((unsigned char) tolower ((char) j));
+ emitPairs[i][j] = emitPairs[i][lj] = emitPairs[li][j] = emitPairs[li][lj] =
+ emitPairs[j][i] = emitPairs[j][li] = emitPairs[lj][i] = emitPairs[lj][li] = exp(pairCounts[j][i]) / totalPairCounts;
+ }
+ emitSingle[i] = emitSingle[li] = exp(singleCounts[i]) / totalSingleCounts;
+ }
+ }
+ }
+
+ /////////////////////////////////////////////////////////////////
+ // ProbabilisticModel::ComputeAlignment()
+ //
+ // Computes an alignment based on the given posterior matrix.
+ // This is done by finding the maximum summing path (or
+ // maximum weight trace) through the posterior matrix. The
+ // final alignment is returned as a pair consisting of:
+ // (1) a string (e.g., XXXBBXXXBBBBBBYYYYBBB) where X's and
+ // denote insertions in one of the two sequences and
+ // B's denote that both sequences are present (i.e.
+ // matches).
+ // (2) a float indicating the sum achieved
+ /////////////////////////////////////////////////////////////////
+
+ pair<SafeVector<char> *, float> ComputeAlignment (int seq1Length, int seq2Length,
+ const VF &posterior) const {
+
+ float *twoRows = new float[(seq2Length+1)*2]; assert (twoRows);
+ float *oldRow = twoRows;
+ float *newRow = twoRows + seq2Length + 1;
+
+ char *tracebackMatrix = new char[(seq1Length+1)*(seq2Length+1)]; assert (tracebackMatrix);
+ char *tracebackPtr = tracebackMatrix;
+
+ VF::const_iterator posteriorPtr = posterior.begin() + seq2Length + 1;
+
+ // initialization
+ for (int i = 0; i <= seq2Length; i++){
+ oldRow[i] = 0;
+ *(tracebackPtr++) = 'L';
+ }
+
+ // fill in matrix
+ for (int i = 1; i <= seq1Length; i++){
+
+ // initialize left column
+ newRow[0] = 0;
+ posteriorPtr++;
+ *(tracebackPtr++) = 'U';
+
+ // fill in rest of row
+ for (int j = 1; j <= seq2Length; j++){
+ ChooseBestOfThree (*(posteriorPtr++) + oldRow[j-1], newRow[j-1], oldRow[j],
+ 'D', 'L', 'U', &newRow[j], tracebackPtr++);
+ }
+
+ // swap rows
+ float *temp = oldRow;
+ oldRow = newRow;
+ newRow = temp;
+ }
+
+ // store best score
+ float total = oldRow[seq2Length];
+ delete [] twoRows;
+
+ // compute traceback
+ SafeVector<char> *alignment = new SafeVector<char>; assert (alignment);
+ int r = seq1Length, c = seq2Length;
+ while (r != 0 || c != 0){
+ char ch = tracebackMatrix[r*(seq2Length+1) + c];
+ switch (ch){
+ case 'L': c--; alignment->push_back ('Y'); break;
+ case 'U': r--; alignment->push_back ('X'); break;
+ case 'D': c--; r--; alignment->push_back ('B'); break;
+ default: assert (false);
+ }
+ }
+
+ delete [] tracebackMatrix;
+
+ reverse (alignment->begin(), alignment->end());
+
+ return make_pair(alignment, total);
+ }
+
+ /////////////////////////////////////////////////////////////////
+ // ProbabilisticModel::ComputeAlignmentWithGapPenalties()
+ //
+ // Similar to ComputeAlignment() except with gap penalties.
+ /////////////////////////////////////////////////////////////////
+
+ pair<SafeVector<char> *, float> ComputeAlignmentWithGapPenalties (MultiSequence *align1,
+ MultiSequence *align2,
+ const VF &posterior, int numSeqs1,
+ int numSeqs2,
+ float gapOpenPenalty,
+ float gapContinuePenalty) const {
+ int seq1Length = align1->GetSequence(0)->GetLength();
+ int seq2Length = align2->GetSequence(0)->GetLength();
+ SafeVector<SafeVector<char>::iterator > dataPtrs1 (align1->GetNumSequences());
+ SafeVector<SafeVector<char>::iterator > dataPtrs2 (align2->GetNumSequences());
+
+ // grab character data
+ for (int i = 0; i < align1->GetNumSequences(); i++)
+ dataPtrs1[i] = align1->GetSequence(i)->GetDataPtr();
+ for (int i = 0; i < align2->GetNumSequences(); i++)
+ dataPtrs2[i] = align2->GetSequence(i)->GetDataPtr();
+
+ // the number of active sequences at any given column is defined to be the
+ // number of non-gap characters in that column; the number of gap opens at
+ // any given column is defined to be the number of gap characters in that
+ // column where the previous character in the respective sequence was not
+ // a gap
+ SafeVector<int> numActive1 (seq1Length+1), numGapOpens1 (seq1Length+1);
+ SafeVector<int> numActive2 (seq2Length+1), numGapOpens2 (seq2Length+1);
+
+ // compute number of active sequences and gap opens for each group
+ for (int i = 0; i < align1->GetNumSequences(); i++){
+ SafeVector<char>::iterator dataPtr = align1->GetSequence(i)->GetDataPtr();
+ numActive1[0] = numGapOpens1[0] = 0;
+ for (int j = 1; j <= seq1Length; j++){
+ if (dataPtr[j] != '-'){
+ numActive1[j]++;
+ numGapOpens1[j] += (j != 1 && dataPtr[j-1] != '-');
+ }
+ }
+ }
+ for (int i = 0; i < align2->GetNumSequences(); i++){
+ SafeVector<char>::iterator dataPtr = align2->GetSequence(i)->GetDataPtr();
+ numActive2[0] = numGapOpens2[0] = 0;
+ for (int j = 1; j <= seq2Length; j++){
+ if (dataPtr[j] != '-'){
+ numActive2[j]++;
+ numGapOpens2[j] += (j != 1 && dataPtr[j-1] != '-');
+ }
+ }
+ }
+
+ VVF openingPenalty1 (numSeqs1+1, VF (numSeqs2+1));
+ VF continuingPenalty1 (numSeqs1+1);
+ VVF openingPenalty2 (numSeqs1+1, VF (numSeqs2+1));
+ VF continuingPenalty2 (numSeqs2+1);
+
+ // precompute penalties
+ for (int i = 0; i <= numSeqs1; i++)
+ for (int j = 0; j <= numSeqs2; j++)
+ openingPenalty1[i][j] = i * (gapOpenPenalty * j + gapContinuePenalty * (numSeqs2 - j));
+ for (int i = 0; i <= numSeqs1; i++)
+ continuingPenalty1[i] = i * gapContinuePenalty * numSeqs2;
+ for (int i = 0; i <= numSeqs2; i++)
+ for (int j = 0; j <= numSeqs1; j++)
+ openingPenalty2[i][j] = i * (gapOpenPenalty * j + gapContinuePenalty * (numSeqs1 - j));
+ for (int i = 0; i <= numSeqs2; i++)
+ continuingPenalty2[i] = i * gapContinuePenalty * numSeqs1;
+
+ float *twoRows = new float[6*(seq2Length+1)]; assert (twoRows);
+ float *oldRowMatch = twoRows;
+ float *newRowMatch = twoRows + (seq2Length+1);
+ float *oldRowInsertX = twoRows + 2*(seq2Length+1);
+ float *newRowInsertX = twoRows + 3*(seq2Length+1);
+ float *oldRowInsertY = twoRows + 4*(seq2Length+1);
+ float *newRowInsertY = twoRows + 5*(seq2Length+1);
+
+ char *tracebackMatrix = new char[3*(seq1Length+1)*(seq2Length+1)]; assert (tracebackMatrix);
+ char *tracebackPtr = tracebackMatrix;
+
+ VF::const_iterator posteriorPtr = posterior.begin() + seq2Length + 1;
+
+ // initialization
+ for (int i = 0; i <= seq2Length; i++){
+ oldRowMatch[i] = oldRowInsertX[i] = (i == 0) ? 0 : LOG_ZERO;
+ oldRowInsertY[i] = (i == 0) ? 0 : oldRowInsertY[i-1] + continuingPenalty2[numActive2[i]];
+ *(tracebackPtr) = *(tracebackPtr+1) = *(tracebackPtr+2) = 'Y';
+ tracebackPtr += 3;
+ }
+
+ // fill in matrix
+ for (int i = 1; i <= seq1Length; i++){
+
+ // initialize left column
+ newRowMatch[0] = newRowInsertY[0] = LOG_ZERO;
+ newRowInsertX[0] = oldRowInsertX[0] + continuingPenalty1[numActive1[i]];
+ posteriorPtr++;
+ *(tracebackPtr) = *(tracebackPtr+1) = *(tracebackPtr+2) = 'X';
+ tracebackPtr += 3;
+
+ // fill in rest of row
+ for (int j = 1; j <= seq2Length; j++){
+
+ // going to MATCH state
+ ChooseBestOfThree (oldRowMatch[j-1],
+ oldRowInsertX[j-1],
+ oldRowInsertY[j-1],
+ 'M', 'X', 'Y', &newRowMatch[j], tracebackPtr++);
+ newRowMatch[j] += *(posteriorPtr++);
+
+ // going to INSERT X state
+ ChooseBestOfThree (oldRowMatch[j] + openingPenalty1[numActive1[i]][numGapOpens2[j]],
+ oldRowInsertX[j] + continuingPenalty1[numActive1[i]],
+ oldRowInsertY[j] + openingPenalty1[numActive1[i]][numGapOpens2[j]],
+ 'M', 'X', 'Y', &newRowInsertX[j], tracebackPtr++);
+
+ // going to INSERT Y state
+ ChooseBestOfThree (newRowMatch[j-1] + openingPenalty2[numActive2[j]][numGapOpens1[i]],
+ newRowInsertX[j-1] + openingPenalty2[numActive2[j]][numGapOpens1[i]],
+ newRowInsertY[j-1] + continuingPenalty2[numActive2[j]],
+ 'M', 'X', 'Y', &newRowInsertY[j], tracebackPtr++);
+ }
+
+ // swap rows
+ float *temp;
+ temp = oldRowMatch; oldRowMatch = newRowMatch; newRowMatch = temp;
+ temp = oldRowInsertX; oldRowInsertX = newRowInsertX; newRowInsertX = temp;
+ temp = oldRowInsertY; oldRowInsertY = newRowInsertY; newRowInsertY = temp;
+ }
+
+ // store best score
+ float total;
+ char matrix;
+ ChooseBestOfThree (oldRowMatch[seq2Length], oldRowInsertX[seq2Length], oldRowInsertY[seq2Length],
+ 'M', 'X', 'Y', &total, &matrix);
+
+ delete [] twoRows;
+
+ // compute traceback
+ SafeVector<char> *alignment = new SafeVector<char>; assert (alignment);
+ int r = seq1Length, c = seq2Length;
+ while (r != 0 || c != 0){
+
+ int offset = (matrix == 'M') ? 0 : (matrix == 'X') ? 1 : 2;
+ char ch = tracebackMatrix[(r*(seq2Length+1) + c) * 3 + offset];
+ switch (matrix){
+ case 'Y': c--; alignment->push_back ('Y'); break;
+ case 'X': r--; alignment->push_back ('X'); break;
+ case 'M': c--; r--; alignment->push_back ('B'); break;
+ default: assert (false);
+ }
+ matrix = ch;
+ }
+
+ delete [] tracebackMatrix;
+
+ reverse (alignment->begin(), alignment->end());
+
+ return make_pair(alignment, 1.0f);
+ }
+
+ /////////////////////////////////////////////////////////////////
+ // ProbabilisticModel::ComputeViterbiAlignment()
+ //
+ // Computes the highest probability pairwise alignment using the
+ // probabilistic model. The final alignment is returned as a
+ // pair consisting of:
+ // (1) a string (e.g., XXXBBXXXBBBBBBYYYYBBB) where X's and
+ // denote insertions in one of the two sequences and
+ // B's denote that both sequences are present (i.e.
+ // matches).
+ // (2) a float containing the log probability of the best
+ // alignment (not used)
+ /////////////////////////////////////////////////////////////////
+
+ pair<SafeVector<char> *, float> ComputeViterbiAlignment (Sequence *seq1, Sequence *seq2) const {
+
+ assert (seq1);
+ assert (seq2);
+
+ const int seq1Length = seq1->GetLength();
+ const int seq2Length = seq2->GetLength();
+
+ // retrieve the points to the beginning of each sequence
+ SafeVector<char>::iterator iter1 = seq1->GetDataPtr();
+ SafeVector<char>::iterator iter2 = seq2->GetDataPtr();
+
+ // create viterbi matrix
+ VF *viterbiPtr = new VF (NumMatrixTypes * (seq1Length+1) * (seq2Length+1), LOG_ZERO);
+ assert (viterbiPtr);
+ VF &viterbi = *viterbiPtr;
+
+ // create traceback matrix
+ VI *tracebackPtr = new VI (NumMatrixTypes * (seq1Length+1) * (seq2Length+1), -1);
+ assert (tracebackPtr);
+ VI &traceback = *tracebackPtr;
+
+ // initialization condition
+ for (int k = 0; k < NumMatrixTypes; k++)
+ viterbi[k] = initialDistribution[k];
+
+ // remember offset for each index combination
+ int ij = 0;
+ int i1j = -seq2Length - 1;
+ int ij1 = -1;
+ int i1j1 = -seq2Length - 2;
+
+ ij *= NumMatrixTypes;
+ i1j *= NumMatrixTypes;
+ ij1 *= NumMatrixTypes;
+ i1j1 *= NumMatrixTypes;
+
+ // compute viterbi scores
+ for (int i = 0; i <= seq1Length; i++){
+ unsigned char c1 = (i == 0) ? '~' : (unsigned char) iter1[i];
+ for (int j = 0; j <= seq2Length; j++){
+ unsigned char c2 = (j == 0) ? '~' : (unsigned char) iter2[j];
+
+ if (i > 0 && j > 0){
+ for (int k = 0; k < NumMatrixTypes; k++){
+ float newVal = viterbi[k + i1j1] + transProb[k][0] + matchProb[c1][c2];
+ if (viterbi[0 + ij] < newVal){
+ viterbi[0 + ij] = newVal;
+ traceback[0 + ij] = k;
+ }
+ }
+ }
+ if (i > 0){
+ for (int k = 0; k < NumInsertStates; k++){
+ float valFromMatch = insProb[c1][k] + viterbi[0 + i1j] + transProb[0][2*k+1];
+ float valFromIns = insProb[c1][k] + viterbi[2*k+1 + i1j] + transProb[2*k+1][2*k+1];
+ if (valFromMatch >= valFromIns){
+ viterbi[2*k+1 + ij] = valFromMatch;
+ traceback[2*k+1 + ij] = 0;
+ }
+ else {
+ viterbi[2*k+1 + ij] = valFromIns;
+ traceback[2*k+1 + ij] = 2*k+1;
+ }
+ }
+ }
+ if (j > 0){
+ for (int k = 0; k < NumInsertStates; k++){
+ float valFromMatch = insProb[c2][k] + viterbi[0 + ij1] + transProb[0][2*k+2];
+ float valFromIns = insProb[c2][k] + viterbi[2*k+2 + ij1] + transProb[2*k+2][2*k+2];
+ if (valFromMatch >= valFromIns){
+ viterbi[2*k+2 + ij] = valFromMatch;
+ traceback[2*k+2 + ij] = 0;
+ }
+ else {
+ viterbi[2*k+2 + ij] = valFromIns;
+ traceback[2*k+2 + ij] = 2*k+2;
+ }
+ }
+ }
+
+ ij += NumMatrixTypes;
+ i1j += NumMatrixTypes;
+ ij1 += NumMatrixTypes;
+ i1j1 += NumMatrixTypes;
+ }
+ }
+
+ // figure out best terminating cell
+ float bestProb = LOG_ZERO;
+ int state = -1;
+ for (int k = 0; k < NumMatrixTypes; k++){
+ float thisProb = viterbi[k + NumMatrixTypes * ((seq1Length+1)*(seq2Length+1) - 1)] + initialDistribution[k];
+ if (bestProb < thisProb){
+ bestProb = thisProb;
+ state = k;
+ }
+ }
+ assert (state != -1);
+
+ delete viterbiPtr;
+
+ // compute traceback
+ SafeVector<char> *alignment = new SafeVector<char>; assert (alignment);
+ int r = seq1Length, c = seq2Length;
+ while (r != 0 || c != 0){
+ int newState = traceback[state + NumMatrixTypes * (r * (seq2Length+1) + c)];
+
+ if (state == 0){ c--; r--; alignment->push_back ('B'); }
+ else if (state % 2 == 1){ r--; alignment->push_back ('X'); }
+ else { c--; alignment->push_back ('Y'); }
+
+ state = newState;
+ }
+
+ delete tracebackPtr;
+
+ reverse (alignment->begin(), alignment->end());
+
+ return make_pair(alignment, bestProb);
+ }
+
+ /////////////////////////////////////////////////////////////////
+ // ProbabilisticModel::BuildPosterior()
+ //
+ // Builds a posterior probability matrix needed to align a pair
+ // of alignments. Mathematically, the returned matrix M is
+ // defined as follows:
+ // M[i,j] = sum sum f(s,t,i,j)
+ // s in align1 t in align2
+ // where
+ // [ P(s[i'] <--> t[j'])
+ // [ if s[i'] is a letter in the ith column of align1 and
+ // [ t[j'] it a letter in the jth column of align2
+ // f(s,t,i,j) = [
+ // [ 0 otherwise
+ //
+ /////////////////////////////////////////////////////////////////
+
+ VF *BuildPosterior (MultiSequence *align1, MultiSequence *align2,
+ const SafeVector<SafeVector<SparseMatrix *> > &sparseMatrices,
+ float cutoff = 0.0f) const {
+ const int seq1Length = align1->GetSequence(0)->GetLength();
+ const int seq2Length = align2->GetSequence(0)->GetLength();
+
+ VF *posteriorPtr = new VF((seq1Length+1) * (seq2Length+1), 0); assert (posteriorPtr);
+ VF &posterior = *posteriorPtr;
+ VF::iterator postPtr = posterior.begin();
+
+ // for each s in align1
+ for (int i = 0; i < align1->GetNumSequences(); i++){
+ int first = align1->GetSequence(i)->GetLabel();
+ SafeVector<int> *mapping1 = align1->GetSequence(i)->GetMapping();
+
+ // for each t in align2
+ for (int j = 0; j < align2->GetNumSequences(); j++){
+ int second = align2->GetSequence(j)->GetLabel();
+ SafeVector<int> *mapping2 = align2->GetSequence(j)->GetMapping();
+
+ if (first < second){
+
+ // get the associated sparse matrix
+ SparseMatrix *matrix = sparseMatrices[first][second];
+
+ for (int ii = 1; ii <= matrix->GetSeq1Length(); ii++){
+ SafeVector<PIF>::iterator row = matrix->GetRowPtr(ii);
+ int base = (*mapping1)[ii] * (seq2Length+1);
+ int rowSize = matrix->GetRowSize(ii);
+
+ // add in all relevant values
+ for (int jj = 0; jj < rowSize; jj++)
+ posterior[base + (*mapping2)[row[jj].first]] += row[jj].second;
+
+ // subtract cutoff
+ for (int jj = 0; jj < matrix->GetSeq2Length(); jj++)
+ posterior[base + (*mapping2)[jj]] -= cutoff;
+ }
+
+ } else {
+
+ // get the associated sparse matrix
+ SparseMatrix *matrix = sparseMatrices[second][first];
+
+ for (int jj = 1; jj <= matrix->GetSeq1Length(); jj++){
+ SafeVector<PIF>::iterator row = matrix->GetRowPtr(jj);
+ int base = (*mapping2)[jj];
+ int rowSize = matrix->GetRowSize(jj);
+
+ // add in all relevant values
+ for (int ii = 0; ii < rowSize; ii++)
+ posterior[base + (*mapping1)[row[ii].first] * (seq2Length + 1)] += row[ii].second;
+
+ // subtract cutoff
+ for (int ii = 0; ii < matrix->GetSeq2Length(); ii++)
+ posterior[base + (*mapping1)[ii] * (seq2Length + 1)] -= cutoff;
+ }
+
+ }
+
+
+ delete mapping2;
+ }
+
+ delete mapping1;
+ }
+
+ return posteriorPtr;
+ }
+};
+
+#endif
git://source.jalview.org / jabaws.git / blobdiff