[med-svn] r270 - trunk/packages/probcons/trunk/debian
Charles Plessy
charles-guest at alioth.debian.org
Wed Apr 18 13:39:19 UTC 2007
Author: charles-guest
Date: 2007-04-18 13:39:18 +0000 (Wed, 18 Apr 2007)
New Revision: 270
Removed:
trunk/packages/probcons/trunk/debian/Main-RNA.cc
Modified:
trunk/packages/probcons/trunk/debian/README.Debian
trunk/packages/probcons/trunk/debian/changelog
trunk/packages/probcons/trunk/debian/rules
Log:
Oops, probcons-RNA was built from obsolete files
Deleted: trunk/packages/probcons/trunk/debian/Main-RNA.cc
===================================================================
--- trunk/packages/probcons/trunk/debian/Main-RNA.cc 2007-04-12 10:59:18 UTC (rev 269)
+++ trunk/packages/probcons/trunk/debian/Main-RNA.cc 2007-04-18 13:39:18 UTC (rev 270)
@@ -1,1445 +0,0 @@
-/////////////////////////////////////////////////////////////////
-// Main.cc
-//
-// Main routines for PROBCONS program.
-/////////////////////////////////////////////////////////////////
-
-#include "SafeVector.h"
-#include "MultiSequence.h"
-#include "Defaults-RNA.h"
-#include "ScoreType.h"
-#include "ProbabilisticModel.h"
-#include "EvolutionaryTree.h"
-#include "SparseMatrix.h"
-#include <string>
-#include <sstream>
-#include <iomanip>
-#include <iostream>
-#include <list>
-#include <set>
-#include <algorithm>
-#include <cstdio>
-#include <cstdlib>
-#include <cerrno>
-#include <iomanip>
-
-string parametersInputFilename = "";
-string parametersOutputFilename = "no training";
-string annotationFilename = "";
-
-bool enableTraining = false;
-bool enableVerbose = false;
-bool enableAllPairs = false;
-bool enableAnnotation = false;
-bool enableViterbi = false;
-bool enableClustalWOutput = false;
-bool enableTrainEmissions = false;
-bool enableAlignOrder = false;
-int numConsistencyReps = 2;
-int numPreTrainingReps = 0;
-int numIterativeRefinementReps = 100;
-
-float cutoff = 0;
-float gapOpenPenalty = 0;
-float gapContinuePenalty = 0;
-
-VF initDistrib (NumMatrixTypes);
-VF gapOpen (2*NumInsertStates);
-VF gapExtend (2*NumInsertStates);
-VVF emitPairs (256, VF (256, 1e-10));
-VF emitSingle (256, 1e-5);
-
-string alphabet = alphabetDefault;
-
-const int MIN_PRETRAINING_REPS = 0;
-const int MAX_PRETRAINING_REPS = 20;
-const int MIN_CONSISTENCY_REPS = 0;
-const int MAX_CONSISTENCY_REPS = 5;
-const int MIN_ITERATIVE_REFINEMENT_REPS = 0;
-const int MAX_ITERATIVE_REFINEMENT_REPS = 1000;
-
-/////////////////////////////////////////////////////////////////
-// Function prototypes
-/////////////////////////////////////////////////////////////////
-
-void PrintHeading();
-void PrintParameters (const char *message, const VF &initDistrib, const VF &gapOpen,
- const VF &gapExtend, const VVF &emitPairs, const VF &emitSingle, const char *filename);
-MultiSequence *DoAlign (MultiSequence *sequence, const ProbabilisticModel &model, VF &initDistrib, VF &gapOpen, VF &gapExtend,
- VVF &emitPairs, VF &emitSingle);
-SafeVector<string> ParseParams (int argc, char **argv);
-void ReadParameters ();
-MultiSequence *ComputeFinalAlignment (const TreeNode *tree, MultiSequence *sequences,
- const SafeVector<SafeVector<SparseMatrix *> > &sparseMatrices,
- const ProbabilisticModel &model);
-MultiSequence *AlignAlignments (MultiSequence *align1, MultiSequence *align2,
- const SafeVector<SafeVector<SparseMatrix *> > &sparseMatrices,
- const ProbabilisticModel &model);
-SafeVector<SafeVector<SparseMatrix *> > DoRelaxation (MultiSequence *sequences,
- SafeVector<SafeVector<SparseMatrix *> > &sparseMatrices);
-void Relax (SparseMatrix *matXZ, SparseMatrix *matZY, VF &posterior);
-void Relax1 (SparseMatrix *matXZ, SparseMatrix *matZY, VF &posterior);
-
-set<int> GetSubtree (const TreeNode *tree);
-void TreeBasedBiPartitioning (const SafeVector<SafeVector<SparseMatrix *> > &sparseMatrices,
- const ProbabilisticModel &model, MultiSequence* &alignment,
- const TreeNode *tree);
-void DoIterativeRefinement (const SafeVector<SafeVector<SparseMatrix *> > &sparseMatrices,
- const ProbabilisticModel &model, MultiSequence* &alignment);
-void WriteAnnotation (MultiSequence *alignment,
- const SafeVector<SafeVector<SparseMatrix *> > &sparseMatrices);
-int ComputeScore (const SafeVector<pair<int, int> > &active,
- const SafeVector<SafeVector<SparseMatrix *> > &sparseMatrices);
-
-
-/////////////////////////////////////////////////////////////////
-// main()
-//
-// Calls all initialization routines and runs the PROBCONS
-// aligner.
-/////////////////////////////////////////////////////////////////
-
-int main (int argc, char **argv){
-
- // print PROBCONS heading
- PrintHeading();
-
- // parse program parameters
- SafeVector<string> sequenceNames = ParseParams (argc, argv);
- ReadParameters();
- PrintParameters ("Using parameter set:", initDistrib, gapOpen, gapExtend, emitPairs, emitSingle, NULL);
-
- // now, we'll process all the files given as input. If we are given
- // several filenames as input, then we'll load all of those sequences
- // simultaneously, as long as we're not training. On the other hand,
- // if we are training, then we'll treat each file as a separate
- // training instance
-
- // if we are training
- if (enableTraining){
-
- // build new model for aligning
- ProbabilisticModel model (initDistrib, gapOpen, gapExtend, emitPairs, emitSingle);
-
- // prepare to average parameters
- for (int i = 0; i < (int) initDistrib.size(); i++) initDistrib[i] = 0;
- for (int i = 0; i < (int) gapOpen.size(); i++) gapOpen[i] = 0;
- for (int i = 0; i < (int) gapExtend.size(); i++) gapExtend[i] = 0;
- if (enableTrainEmissions){
- for (int i = 0; i < (int) emitPairs.size(); i++)
- for (int j = 0; j < (int) emitPairs[i].size(); j++) emitPairs[i][j] = 0;
- for (int i = 0; i < (int) emitSingle.size(); i++) emitSingle[i] = 0;
- }
-
- // align each file individually
- for (int i = 0; i < (int) sequenceNames.size(); i++){
-
- VF thisInitDistrib (NumMatrixTypes);
- VF thisGapOpen (2*NumInsertStates);
- VF thisGapExtend (2*NumInsertStates);
- VVF thisEmitPairs (256, VF (256, 1e-10));
- VF thisEmitSingle (256, 1e-5);
-
- // load sequence file
- MultiSequence *sequences = new MultiSequence(); assert (sequences);
- cerr << "Loading sequence file: " << sequenceNames[i] << endl;
- sequences->LoadMFA (sequenceNames[i], true);
-
- // align sequences
- DoAlign (sequences, model, thisInitDistrib, thisGapOpen, thisGapExtend, thisEmitPairs, thisEmitSingle);
-
- // add in contribution of the derived parameters
- for (int i = 0; i < (int) initDistrib.size(); i++) initDistrib[i] += thisInitDistrib[i];
- for (int i = 0; i < (int) gapOpen.size(); i++) gapOpen[i] += thisGapOpen[i];
- for (int i = 0; i < (int) gapExtend.size(); i++) gapExtend[i] += thisGapExtend[i];
- if (enableTrainEmissions){
- for (int i = 0; i < (int) emitPairs.size(); i++)
- for (int j = 0; j < (int) emitPairs[i].size(); j++) emitPairs[i][j] += thisEmitPairs[i][j];
- for (int i = 0; i < (int) emitSingle.size(); i++) emitSingle[i] += thisEmitSingle[i];
- }
-
- delete sequences;
- }
-
- // compute new parameters and print them out
- for (int i = 0; i < (int) initDistrib.size(); i++) initDistrib[i] /= (int) sequenceNames.size();
- for (int i = 0; i < (int) gapOpen.size(); i++) gapOpen[i] /= (int) sequenceNames.size();
- for (int i = 0; i < (int) gapExtend.size(); i++) gapExtend[i] /= (int) sequenceNames.size();
- if (enableTrainEmissions){
- for (int i = 0; i < (int) emitPairs.size(); i++)
- for (int j = 0; j < (int) emitPairs[i].size(); j++) emitPairs[i][j] /= (int) sequenceNames.size();
- for (int i = 0; i < (int) emitSingle.size(); i++) emitSingle[i] /= sequenceNames.size();
- }
-
- PrintParameters ("Trained parameter set:",
- initDistrib, gapOpen, gapExtend, emitPairs, emitSingle,
- parametersOutputFilename.c_str());
- }
-
- // if we are not training, we must simply want to align some sequences
- else {
-
- // load all files together
- MultiSequence *sequences = new MultiSequence(); assert (sequences);
- for (int i = 0; i < (int) sequenceNames.size(); i++){
- cerr << "Loading sequence file: " << sequenceNames[i] << endl;
- sequences->LoadMFA (sequenceNames[i], true);
- }
-
- // do all "pre-training" repetitions first
- for (int ct = 0; ct < numPreTrainingReps; ct++){
- enableTraining = true;
-
- // build new model for aligning
- ProbabilisticModel model (initDistrib, gapOpen, gapExtend,
- emitPairs, emitSingle);
-
- // do initial alignments
- DoAlign (sequences, model, initDistrib, gapOpen, gapExtend, emitPairs, emitSingle);
-
- // print new parameters
- PrintParameters ("Recomputed parameter set:", initDistrib, gapOpen, gapExtend, emitPairs, emitSingle, NULL);
-
- enableTraining = false;
- }
-
- // now, we can perform the alignments and write them out
- MultiSequence *alignment = DoAlign (sequences,
- ProbabilisticModel (initDistrib, gapOpen, gapExtend, emitPairs, emitSingle),
- initDistrib, gapOpen, gapExtend, emitPairs, emitSingle);
-
- if (!enableAllPairs){
- if (enableClustalWOutput)
- alignment->WriteALN (cout);
- else
- alignment->WriteMFA (cout);
- }
- delete alignment;
- delete sequences;
-
- }
-}
-
-/////////////////////////////////////////////////////////////////
-// PrintHeading()
-//
-// Prints heading for PROBCONS program.
-/////////////////////////////////////////////////////////////////
-
-void PrintHeading (){
- cerr << endl
- << "PROBCONS version " << VERSION << " - align multiple protein sequences and print to standard output" << endl
- << "Written by Chuong Do" << endl
- << endl;
-}
-
-/////////////////////////////////////////////////////////////////
-// PrintParameters()
-//
-// Prints PROBCONS parameters to STDERR. If a filename is
-// specified, then the parameters are also written to the file.
-/////////////////////////////////////////////////////////////////
-
-void PrintParameters (const char *message, const VF &initDistrib, const VF &gapOpen,
- const VF &gapExtend, const VVF &emitPairs, const VF &emitSingle, const char *filename){
-
- // print parameters to the screen
- cerr << message << endl
- << " initDistrib[] = { ";
- for (int i = 0; i < NumMatrixTypes; i++) cerr << setprecision (10) << initDistrib[i] << " ";
- cerr << "}" << endl
- << " gapOpen[] = { ";
- for (int i = 0; i < NumInsertStates*2; i++) cerr << setprecision (10) << gapOpen[i] << " ";
- cerr << "}" << endl
- << " gapExtend[] = { ";
- for (int i = 0; i < NumInsertStates*2; i++) cerr << setprecision (10) << gapExtend[i] << " ";
- cerr << "}" << endl
- << endl;
-
- /*
- for (int i = 0; i < 5; i++){
- for (int j = 0; j <= i; j++){
- cerr << emitPairs[(unsigned char) alphabet[i]][(unsigned char) alphabet[j]] << " ";
- }
- cerr << endl;
- }*/
-
- // if a file name is specified
- if (filename){
-
- // attempt to open the file for writing
- FILE *file = fopen (filename, "w");
- if (!file){
- cerr << "ERROR: Unable to write parameter file: " << filename << endl;
- exit (1);
- }
-
- // if successful, then write the parameters to the file
- for (int i = 0; i < NumMatrixTypes; i++) fprintf (file, "%.10f ", initDistrib[i]); fprintf (file, "\n");
- for (int i = 0; i < 2*NumInsertStates; i++) fprintf (file, "%.10f ", gapOpen[i]); fprintf (file, "\n");
- for (int i = 0; i < 2*NumInsertStates; i++) fprintf (file, "%.10f ", gapExtend[i]); fprintf (file, "\n");
- fprintf (file, "%s\n", alphabet.c_str());
- for (int i = 0; i < (int) alphabet.size(); i++){
- for (int j = 0; j <= i; j++)
- fprintf (file, "%.10f ", emitPairs[(unsigned char) alphabet[i]][(unsigned char) alphabet[j]]);
- fprintf (file, "\n");
- }
- for (int i = 0; i < (int) alphabet.size(); i++)
- fprintf (file, "%.10f ", emitSingle[(unsigned char) alphabet[i]]);
- fprintf (file, "\n");
- fclose (file);
- }
-}
-
-/////////////////////////////////////////////////////////////////
-// DoAlign()
-//
-// First computes all pairwise posterior probability matrices.
-// Then, computes new parameters if training, or a final
-// alignment, otherwise.
-/////////////////////////////////////////////////////////////////
-
-MultiSequence *DoAlign (MultiSequence *sequences, const ProbabilisticModel &model, VF &initDistrib, VF &gapOpen,
- VF &gapExtend, VVF &emitPairs, VF &emitSingle){
-
- assert (sequences);
-
- const int numSeqs = sequences->GetNumSequences();
- VVF distances (numSeqs, VF (numSeqs, 0));
- SafeVector<SafeVector<SparseMatrix *> > sparseMatrices (numSeqs, SafeVector<SparseMatrix *>(numSeqs, NULL));
-
-
-
- if (enableTraining){
- // prepare to average parameters
- for (int i = 0; i < (int) initDistrib.size(); i++) initDistrib[i] = 0;
- for (int i = 0; i < (int) gapOpen.size(); i++) gapOpen[i] = 0;
- for (int i = 0; i < (int) gapExtend.size(); i++) gapExtend[i] = 0;
- if (enableTrainEmissions){
- for (int i = 0; i < (int) emitPairs.size(); i++)
- for (int j = 0; j < (int) emitPairs[i].size(); j++) emitPairs[i][j] = 0;
- for (int i = 0; i < (int) emitSingle.size(); i++) emitSingle[i] = 0;
- }
- }
-
- // skip posterior calculations if we just want to do Viterbi alignments
- if (!enableViterbi){
-
- // do all pairwise alignments for posterior probability matrices
- for (int a = 0; a < numSeqs-1; a++){
- for (int b = a+1; b < numSeqs; b++){
- Sequence *seq1 = sequences->GetSequence (a);
- Sequence *seq2 = sequences->GetSequence (b);
-
- // verbose output
- if (enableVerbose)
- cerr << "Computing posterior matrix: (" << a+1 << ") " << seq1->GetHeader() << " vs. "
- << "(" << b+1 << ") " << seq2->GetHeader() << " -- ";
-
- // compute forward and backward probabilities
- VF *forward = model.ComputeForwardMatrix (seq1, seq2); assert (forward);
- VF *backward = model.ComputeBackwardMatrix (seq1, seq2); assert (backward);
-
- // if we are training, then we'll simply want to compute the
- // expected counts for each region within the matrix separately;
- // otherwise, we'll need to put all of the regions together and
- // assemble a posterior probability match matrix
-
- // so, if we're training
- if (enableTraining){
-
- // compute new parameters
- VF thisInitDistrib (NumMatrixTypes);
- VF thisGapOpen (2*NumInsertStates);
- VF thisGapExtend (2*NumInsertStates);
- VVF thisEmitPairs (256, VF (256, 1e-10));
- VF thisEmitSingle (256, 1e-5);
-
- model.ComputeNewParameters (seq1, seq2, *forward, *backward, thisInitDistrib, thisGapOpen, thisGapExtend, thisEmitPairs, thisEmitSingle, enableTrainEmissions);
-
- // add in contribution of the derived parameters
- for (int i = 0; i < (int) initDistrib.size(); i++) initDistrib[i] += thisInitDistrib[i];
- for (int i = 0; i < (int) gapOpen.size(); i++) gapOpen[i] += thisGapOpen[i];
- for (int i = 0; i < (int) gapExtend.size(); i++) gapExtend[i] += thisGapExtend[i];
- if (enableTrainEmissions){
- for (int i = 0; i < (int) emitPairs.size(); i++)
- for (int j = 0; j < (int) emitPairs[i].size(); j++) emitPairs[i][j] += thisEmitPairs[i][j];
- for (int i = 0; i < (int) emitSingle.size(); i++) emitSingle[i] += thisEmitSingle[i];
- }
-
- // let us know that we're done.
- if (enableVerbose) cerr << "done." << endl;
- }
- else {
-
- // compute posterior probability matrix
- VF *posterior = model.ComputePosteriorMatrix (seq1, seq2, *forward, *backward); assert (posterior);
-
- // compute sparse representations
- sparseMatrices[a][b] = new SparseMatrix (seq1->GetLength(), seq2->GetLength(), *posterior);
- sparseMatrices[b][a] = NULL;
-
- if (!enableAllPairs){
- // perform the pairwise sequence alignment
- pair<SafeVector<char> *, float> alignment = model.ComputeAlignment (seq1->GetLength(),
- seq2->GetLength(),
- *posterior);
-
- // compute "expected accuracy" distance for evolutionary tree computation
- float distance = alignment.second / min (seq1->GetLength(), seq2->GetLength());
- distances[a][b] = distances[b][a] = distance;
-
- if (enableVerbose)
- cerr << setprecision (10) << distance << endl;
-
- delete alignment.first;
- }
- else {
- // let us know that we're done.
- if (enableVerbose) cerr << "done." << endl;
- }
-
- delete posterior;
- }
-
- delete forward;
- delete backward;
- }
- }
- }
-
- // now average out parameters derived
- if (enableTraining){
-
- // compute new parameters
- for (int i = 0; i < (int) initDistrib.size(); i++) initDistrib[i] /= numSeqs * (numSeqs - 1) / 2;
- for (int i = 0; i < (int) gapOpen.size(); i++) gapOpen[i] /= numSeqs * (numSeqs - 1) / 2;
- for (int i = 0; i < (int) gapExtend.size(); i++) gapExtend[i] /= numSeqs * (numSeqs - 1) / 2;
-
- if (enableTrainEmissions){
- for (int i = 0; i < (int) emitPairs.size(); i++)
- for (int j = 0; j < (int) emitPairs[i].size(); j++) emitPairs[i][j] /= numSeqs * (numSeqs - 1) / 2;
- for (int i = 0; i < (int) emitSingle.size(); i++) emitSingle[i] /= numSeqs * (numSeqs - 1) / 2;
- }
- }
-
- // see if we still want to do some alignments
- else {
-
- if (!enableViterbi){
-
- // perform the consistency transformation the desired number of times
- for (int r = 0; r < numConsistencyReps; r++){
- SafeVector<SafeVector<SparseMatrix *> > newSparseMatrices = DoRelaxation (sequences, sparseMatrices);
-
- // now replace the old posterior matrices
- for (int i = 0; i < numSeqs; i++){
- for (int j = 0; j < numSeqs; j++){
- delete sparseMatrices[i][j];
- sparseMatrices[i][j] = newSparseMatrices[i][j];
- }
- }
- }
- }
-
- MultiSequence *finalAlignment = NULL;
-
- if (enableAllPairs){
- for (int a = 0; a < numSeqs-1; a++){
- for (int b = a+1; b < numSeqs; b++){
- Sequence *seq1 = sequences->GetSequence (a);
- Sequence *seq2 = sequences->GetSequence (b);
-
- if (enableVerbose)
- cerr << "Performing pairwise alignment: (" << a+1 << ") " << seq1->GetHeader() << " vs. "
- << "(" << b+1 << ") " << seq2->GetHeader() << " -- ";
-
-
- // perform the pairwise sequence alignment
- pair<SafeVector<char> *, float> alignment;
- if (enableViterbi)
- alignment = model.ComputeViterbiAlignment (seq1, seq2);
- else {
-
- // build posterior matrix
- VF *posterior = sparseMatrices[a][b]->GetPosterior(); assert (posterior);
- int length = (seq1->GetLength() + 1) * (seq2->GetLength() + 1);
- for (int i = 0; i < length; i++) (*posterior)[i] -= cutoff;
-
- alignment = model.ComputeAlignment (seq1->GetLength(), seq2->GetLength(), *posterior);
- delete posterior;
- }
-
- // write pairwise alignments
- string name = seq1->GetHeader() + "-" + seq2->GetHeader() + (enableClustalWOutput ? ".aln" : ".fasta");
- ofstream outfile (name.c_str());
-
- MultiSequence *result = new MultiSequence();
- result->AddSequence (seq1->AddGaps(alignment.first, 'X'));
- result->AddSequence (seq2->AddGaps(alignment.first, 'Y'));
- if (enableClustalWOutput)
- result->WriteALN (outfile);
- else
- result->WriteMFA (outfile);
-
- outfile.close();
-
- delete alignment.first;
- }
- }
- }
-
- // now if we still need to do a final multiple alignment
- else {
-
- if (enableVerbose)
- cerr << endl;
-
- // compute the evolutionary tree
- TreeNode *tree = TreeNode::ComputeTree (distances);
-
- tree->Print (cerr, sequences);
- cerr << endl;
-
- // make the final alignment
- finalAlignment = ComputeFinalAlignment (tree, sequences, sparseMatrices, model);
-
- // build annotation
- if (enableAnnotation){
- WriteAnnotation (finalAlignment, sparseMatrices);
- }
-
- delete tree;
- }
-
- if (!enableViterbi){
- // delete sparse matrices
- for (int a = 0; a < numSeqs-1; a++){
- for (int b = a+1; b < numSeqs; b++){
- delete sparseMatrices[a][b];
- delete sparseMatrices[b][a];
- }
- }
- }
-
- return finalAlignment;
- }
-
- return NULL;
-}
-
-/////////////////////////////////////////////////////////////////
-// GetInteger()
-//
-// Attempts to parse an integer from the character string given.
-// Returns true only if no parsing error occurs.
-/////////////////////////////////////////////////////////////////
-
-bool GetInteger (char *data, int *val){
- char *endPtr;
- long int retVal;
-
- assert (val);
-
- errno = 0;
- retVal = strtol (data, &endPtr, 0);
- if (retVal == 0 && (errno != 0 || data == endPtr)) return false;
- if (errno != 0 && (retVal == LONG_MAX || retVal == LONG_MIN)) return false;
- if (retVal < (long) INT_MIN || retVal > (long) INT_MAX) return false;
- *val = (int) retVal;
- return true;
-}
-
-/////////////////////////////////////////////////////////////////
-// GetFloat()
-//
-// Attempts to parse a float from the character string given.
-// Returns true only if no parsing error occurs.
-/////////////////////////////////////////////////////////////////
-
-bool GetFloat (char *data, float *val){
- char *endPtr;
- double retVal;
-
- assert (val);
-
- errno = 0;
- retVal = strtod (data, &endPtr);
- if (retVal == 0 && (errno != 0 || data == endPtr)) return false;
- if (errno != 0 && (retVal >= 1000000.0 || retVal <= -1000000.0)) return false;
- *val = (float) retVal;
- return true;
-}
-
-/////////////////////////////////////////////////////////////////
-// ParseParams()
-//
-// Parse all command-line options.
-/////////////////////////////////////////////////////////////////
-
-SafeVector<string> ParseParams (int argc, char **argv){
-
- if (argc < 2){
-
- cerr << "PROBCONS comes with ABSOLUTELY NO WARRANTY. This is free software, and" << endl
- << "you are welcome to redistribute it under certain conditions. See the" << endl
- << "file COPYING.txt for details." << endl
- << endl
- << "Usage:" << endl
- << " probcons [OPTION]... [MFAFILE]..." << endl
- << endl
- << "Description:" << endl
- << " Align sequences in MFAFILE(s) and print result to standard output" << endl
- << endl
- << " -clustalw" << endl
- << " use CLUSTALW output format instead of MFA" << endl
- << endl
- << " -c, --consistency REPS" << endl
- << " use " << MIN_CONSISTENCY_REPS << " <= REPS <= " << MAX_CONSISTENCY_REPS
- << " (default: " << numConsistencyReps << ") passes of consistency transformation" << endl
- << endl
- << " -ir, --iterative-refinement REPS" << endl
- << " use " << MIN_ITERATIVE_REFINEMENT_REPS << " <= REPS <= " << MAX_ITERATIVE_REFINEMENT_REPS
- << " (default: " << numIterativeRefinementReps << ") passes of iterative-refinement" << endl
- << endl
- << " -pre, --pre-training REPS" << endl
- << " use " << MIN_PRETRAINING_REPS << " <= REPS <= " << MAX_PRETRAINING_REPS
- << " (default: " << numPreTrainingReps << ") rounds of pretraining" << endl
- << endl
- << " -pairs" << endl
- << " generate all-pairs pairwise alignments" << endl
- << endl
- << " -viterbi" << endl
- << " use Viterbi algorithm to generate all pairs (automatically enables -pairs)" << endl
- << endl
- << " -v, --verbose" << endl
- << " report progress while aligning (default: " << (enableVerbose ? "on" : "off") << ")" << endl
- << endl
- << " -annot FILENAME" << endl
- << " write annotation for multiple alignment to FILENAME" << endl
- << endl
- << " -t, --train FILENAME" << endl
- << " compute EM transition probabilities, store in FILENAME (default: "
- << parametersOutputFilename << ")" << endl
- << endl
- << " -e, --emissions" << endl
- << " also reestimate emission probabilities (default: "
- << (enableTrainEmissions ? "on" : "off") << ")" << endl
- << endl
- << " -p, --paramfile FILENAME" << endl
- << " read parameters from FILENAME (default: "
- << parametersInputFilename << ")" << endl
- << endl
- << " -a, --alignment-order" << endl
- << " print sequences in alignment order rather than input order (default: "
- << (enableAlignOrder ? "on" : "off") << ")" << endl
- << endl;
- // << " -go, --gap-open VALUE" << endl
- // << " gap opening penalty of VALUE <= 0 (default: " << gapOpenPenalty << ")" << endl
- // << endl
- // << " -ge, --gap-extension VALUE" << endl
- // << " gap extension penalty of VALUE <= 0 (default: " << gapContinuePenalty << ")" << endl
- // << endl
- // << " -co, --cutoff CUTOFF" << endl
- // << " subtract 0 <= CUTOFF <= 1 (default: " << cutoff << ") from all posterior values before final alignment" << endl
- // << endl
-
- exit (1);
- }
-
- SafeVector<string> sequenceNames;
- int tempInt;
- float tempFloat;
-
- for (int i = 1; i < argc; i++){
- if (argv[i][0] == '-'){
-
- // training
- if (!strcmp (argv[i], "-t") || !strcmp (argv[i], "--train")){
- enableTraining = true;
- if (i < argc - 1)
- parametersOutputFilename = string (argv[++i]);
- else {
- cerr << "ERROR: Filename expected for option " << argv[i] << endl;
- exit (1);
- }
- }
-
- // emission training
- else if (!strcmp (argv[i], "-e") || !strcmp (argv[i], "--emissions")){
- enableTrainEmissions = true;
- }
-
- // parameter file
- else if (!strcmp (argv[i], "-p") || !strcmp (argv[i], "--paramfile")){
- if (i < argc - 1)
- parametersInputFilename = string (argv[++i]);
- else {
- cerr << "ERROR: Filename expected for option " << argv[i] << endl;
- exit (1);
- }
- }
-
- // number of consistency transformations
- else if (!strcmp (argv[i], "-c") || !strcmp (argv[i], "--consistency")){
- if (i < argc - 1){
- if (!GetInteger (argv[++i], &tempInt)){
- cerr << "ERROR: Invalid integer following option " << argv[i-1] << ": " << argv[i] << endl;
- exit (1);
- }
- else {
- if (tempInt < MIN_CONSISTENCY_REPS || tempInt > MAX_CONSISTENCY_REPS){
- cerr << "ERROR: For option " << argv[i-1] << ", integer must be between "
- << MIN_CONSISTENCY_REPS << " and " << MAX_CONSISTENCY_REPS << "." << endl;
- exit (1);
- }
- else
- numConsistencyReps = tempInt;
- }
- }
- else {
- cerr << "ERROR: Integer expected for option " << argv[i] << endl;
- exit (1);
- }
- }
-
- // number of randomized partitioning iterative refinement passes
- else if (!strcmp (argv[i], "-ir") || !strcmp (argv[i], "--iterative-refinement")){
- if (i < argc - 1){
- if (!GetInteger (argv[++i], &tempInt)){
- cerr << "ERROR: Invalid integer following option " << argv[i-1] << ": " << argv[i] << endl;
- exit (1);
- }
- else {
- if (tempInt < MIN_ITERATIVE_REFINEMENT_REPS || tempInt > MAX_ITERATIVE_REFINEMENT_REPS){
- cerr << "ERROR: For option " << argv[i-1] << ", integer must be between "
- << MIN_ITERATIVE_REFINEMENT_REPS << " and " << MAX_ITERATIVE_REFINEMENT_REPS << "." << endl;
- exit (1);
- }
- else
- numIterativeRefinementReps = tempInt;
- }
- }
- else {
- cerr << "ERROR: Integer expected for option " << argv[i] << endl;
- exit (1);
- }
- }
-
- // number of EM pre-training rounds
- else if (!strcmp (argv[i], "-pre") || !strcmp (argv[i], "--pre-training")){
- if (i < argc - 1){
- if (!GetInteger (argv[++i], &tempInt)){
- cerr << "ERROR: Invalid integer following option " << argv[i-1] << ": " << argv[i] << endl;
- exit (1);
- }
- else {
- if (tempInt < MIN_PRETRAINING_REPS || tempInt > MAX_PRETRAINING_REPS){
- cerr << "ERROR: For option " << argv[i-1] << ", integer must be between "
- << MIN_PRETRAINING_REPS << " and " << MAX_PRETRAINING_REPS << "." << endl;
- exit (1);
- }
- else
- numPreTrainingReps = tempInt;
- }
- }
- else {
- cerr << "ERROR: Integer expected for option " << argv[i] << endl;
- exit (1);
- }
- }
-
- // gap open penalty
- else if (!strcmp (argv[i], "-go") || !strcmp (argv[i], "--gap-open")){
- if (i < argc - 1){
- if (!GetFloat (argv[++i], &tempFloat)){
- cerr << "ERROR: Invalid floating-point value following option " << argv[i-1] << ": " << argv[i] << endl;
- exit (1);
- }
- else {
- if (tempFloat > 0){
- cerr << "ERROR: For option " << argv[i-1] << ", floating-point value must not be positive." << endl;
- exit (1);
- }
- else
- gapOpenPenalty = tempFloat;
- }
- }
- else {
- cerr << "ERROR: Floating-point value expected for option " << argv[i] << endl;
- exit (1);
- }
- }
-
- // gap extension penalty
- else if (!strcmp (argv[i], "-ge") || !strcmp (argv[i], "--gap-extension")){
- if (i < argc - 1){
- if (!GetFloat (argv[++i], &tempFloat)){
- cerr << "ERROR: Invalid floating-point value following option " << argv[i-1] << ": " << argv[i] << endl;
- exit (1);
- }
- else {
- if (tempFloat > 0){
- cerr << "ERROR: For option " << argv[i-1] << ", floating-point value must not be positive." << endl;
- exit (1);
- }
- else
- gapContinuePenalty = tempFloat;
- }
- }
- else {
- cerr << "ERROR: Floating-point value expected for option " << argv[i] << endl;
- exit (1);
- }
- }
-
- // all-pairs pairwise alignments
- else if (!strcmp (argv[i], "-pairs")){
- enableAllPairs = true;
- }
-
- // all-pairs pairwise Viterbi alignments
- else if (!strcmp (argv[i], "-viterbi")){
- enableAllPairs = true;
- enableViterbi = true;
- }
-
- // annotation files
- else if (!strcmp (argv[i], "-annot")){
- enableAnnotation = true;
- if (i < argc - 1)
- annotationFilename = argv[++i];
- else {
- cerr << "ERROR: FILENAME expected for option " << argv[i] << endl;
- exit (1);
- }
- }
-
- // clustalw output format
- else if (!strcmp (argv[i], "-clustalw")){
- enableClustalWOutput = true;
- }
-
- // cutoff
- else if (!strcmp (argv[i], "-co") || !strcmp (argv[i], "--cutoff")){
- if (i < argc - 1){
- if (!GetFloat (argv[++i], &tempFloat)){
- cerr << "ERROR: Invalid floating-point value following option " << argv[i-1] << ": " << argv[i] << endl;
- exit (1);
- }
- else {
- if (tempFloat < 0 || tempFloat > 1){
- cerr << "ERROR: For option " << argv[i-1] << ", floating-point value must be between 0 and 1." << endl;
- exit (1);
- }
- else
- cutoff = tempFloat;
- }
- }
- else {
- cerr << "ERROR: Floating-point value expected for option " << argv[i] << endl;
- exit (1);
- }
- }
-
- // verbose reporting
- else if (!strcmp (argv[i], "-v") || !strcmp (argv[i], "--verbose")){
- enableVerbose = true;
- }
-
- // alignment order
- else if (!strcmp (argv[i], "-a") || !strcmp (argv[i], "--alignment-order")){
- enableAlignOrder = true;
- }
-
- // bad arguments
- else {
- cerr << "ERROR: Unrecognized option: " << argv[i] << endl;
- exit (1);
- }
- }
- else {
- sequenceNames.push_back (string (argv[i]));
- }
- }
-
- if (enableTrainEmissions && !enableTraining){
- cerr << "ERROR: Training emissions (-e) requires training (-t)" << endl;
- exit (1);
- }
-
- return sequenceNames;
-}
-
-/////////////////////////////////////////////////////////////////
-// ReadParameters()
-//
-// Read initial distribution, transition, and emission
-// parameters from a file.
-/////////////////////////////////////////////////////////////////
-
-void ReadParameters (){
-
- ifstream data;
-
- emitPairs = VVF (256, VF (256, 1e-10));
- emitSingle = VF (256, 1e-5);
-
- // read initial state distribution and transition parameters
- if (parametersInputFilename == string ("")){
- if (NumInsertStates == 1){
- for (int i = 0; i < NumMatrixTypes; i++) initDistrib[i] = initDistrib1Default[i];
- for (int i = 0; i < 2*NumInsertStates; i++) gapOpen[i] = gapOpen1Default[i];
- for (int i = 0; i < 2*NumInsertStates; i++) gapExtend[i] = gapExtend1Default[i];
- }
- else if (NumInsertStates == 2){
- for (int i = 0; i < NumMatrixTypes; i++) initDistrib[i] = initDistrib2Default[i];
- for (int i = 0; i < 2*NumInsertStates; i++) gapOpen[i] = gapOpen2Default[i];
- for (int i = 0; i < 2*NumInsertStates; i++) gapExtend[i] = gapExtend2Default[i];
- }
- else {
- cerr << "ERROR: No default initial distribution/parameter settings exist" << endl
- << " for " << NumInsertStates << " pairs of insert states. Use --paramfile." << endl;
- exit (1);
- }
-
- alphabet = alphabetDefault;
-
- for (int i = 0; i < (int) alphabet.length(); i++){
- emitSingle[(unsigned char) tolower(alphabet[i])] = emitSingleDefault[i];
- emitSingle[(unsigned char) toupper(alphabet[i])] = emitSingleDefault[i];
- for (int j = 0; j <= i; j++){
- emitPairs[(unsigned char) tolower(alphabet[i])][(unsigned char) tolower(alphabet[j])] = emitPairsDefault[i][j];
- emitPairs[(unsigned char) tolower(alphabet[i])][(unsigned char) toupper(alphabet[j])] = emitPairsDefault[i][j];
- emitPairs[(unsigned char) toupper(alphabet[i])][(unsigned char) tolower(alphabet[j])] = emitPairsDefault[i][j];
- emitPairs[(unsigned char) toupper(alphabet[i])][(unsigned char) toupper(alphabet[j])] = emitPairsDefault[i][j];
- emitPairs[(unsigned char) tolower(alphabet[j])][(unsigned char) tolower(alphabet[i])] = emitPairsDefault[i][j];
- emitPairs[(unsigned char) tolower(alphabet[j])][(unsigned char) toupper(alphabet[i])] = emitPairsDefault[i][j];
- emitPairs[(unsigned char) toupper(alphabet[j])][(unsigned char) tolower(alphabet[i])] = emitPairsDefault[i][j];
- emitPairs[(unsigned char) toupper(alphabet[j])][(unsigned char) toupper(alphabet[i])] = emitPairsDefault[i][j];
- }
- }
- }
- else {
- data.open (parametersInputFilename.c_str());
- if (data.fail()){
- cerr << "ERROR: Unable to read parameter file: " << parametersInputFilename << endl;
- exit (1);
- }
-
- string line[3];
- for (int i = 0; i < 3; i++){
- if (!getline (data, line[i])){
- cerr << "ERROR: Unable to read transition parameters from parameter file: " << parametersInputFilename << endl;
- exit (1);
- }
- }
- istringstream data2;
- data2.clear(); data2.str (line[0]); for (int i = 0; i < NumMatrixTypes; i++) data2 >> initDistrib[i];
- data2.clear(); data2.str (line[1]); for (int i = 0; i < 2*NumInsertStates; i++) data2 >> gapOpen[i];
- data2.clear(); data2.str (line[2]); for (int i = 0; i < 2*NumInsertStates; i++) data2 >> gapExtend[i];
-
- if (!getline (data, line[0])){
- cerr << "ERROR: Unable to read alphabet from scoring matrix file: " << parametersInputFilename << endl;
- exit (1);
- }
-
- // read alphabet as concatenation of all characters on alphabet line
- alphabet = "";
- string token;
- data2.clear(); data2.str (line[0]); while (data2 >> token) alphabet += token;
-
- for (int i = 0; i < (int) alphabet.size(); i++){
- for (int j = 0; j <= i; j++){
- float val;
- data >> val;
- emitPairs[(unsigned char) tolower(alphabet[i])][(unsigned char) tolower(alphabet[j])] = val;
- emitPairs[(unsigned char) tolower(alphabet[i])][(unsigned char) toupper(alphabet[j])] = val;
- emitPairs[(unsigned char) toupper(alphabet[i])][(unsigned char) tolower(alphabet[j])] = val;
- emitPairs[(unsigned char) toupper(alphabet[i])][(unsigned char) toupper(alphabet[j])] = val;
- emitPairs[(unsigned char) tolower(alphabet[j])][(unsigned char) tolower(alphabet[i])] = val;
- emitPairs[(unsigned char) tolower(alphabet[j])][(unsigned char) toupper(alphabet[i])] = val;
- emitPairs[(unsigned char) toupper(alphabet[j])][(unsigned char) tolower(alphabet[i])] = val;
- emitPairs[(unsigned char) toupper(alphabet[j])][(unsigned char) toupper(alphabet[i])] = val;
- }
- }
-
- for (int i = 0; i < (int) alphabet.size(); i++){
- float val;
- data >> val;
- emitSingle[(unsigned char) tolower(alphabet[i])] = val;
- emitSingle[(unsigned char) toupper(alphabet[i])] = val;
- }
- data.close();
- }
-}
-
-/////////////////////////////////////////////////////////////////
-// ProcessTree()
-//
-// Process the tree recursively. Returns the aligned sequences
-// corresponding to a node or leaf of the tree.
-/////////////////////////////////////////////////////////////////
-
-MultiSequence *ProcessTree (const TreeNode *tree, MultiSequence *sequences,
- const SafeVector<SafeVector<SparseMatrix *> > &sparseMatrices,
- const ProbabilisticModel &model){
- MultiSequence *result;
-
- // check if this is a node of the alignment tree
- if (tree->GetSequenceLabel() == -1){
- MultiSequence *alignLeft = ProcessTree (tree->GetLeftChild(), sequences, sparseMatrices, model);
- MultiSequence *alignRight = ProcessTree (tree->GetRightChild(), sequences, sparseMatrices, model);
-
- assert (alignLeft);
- assert (alignRight);
-
- result = AlignAlignments (alignLeft, alignRight, sparseMatrices, model);
- assert (result);
-
- delete alignLeft;
- delete alignRight;
- }
-
- // otherwise, this is a leaf of the alignment tree
- else {
- result = new MultiSequence(); assert (result);
- result->AddSequence (sequences->GetSequence(tree->GetSequenceLabel())->Clone());
- }
-
- return result;
-}
-
-/////////////////////////////////////////////////////////////////
-// ComputeFinalAlignment()
-//
-// Compute the final alignment by calling ProcessTree(), then
-// performing iterative refinement as needed.
-/////////////////////////////////////////////////////////////////
-
-MultiSequence *ComputeFinalAlignment (const TreeNode *tree, MultiSequence *sequences,
- const SafeVector<SafeVector<SparseMatrix *> > &sparseMatrices,
- const ProbabilisticModel &model){
-
- MultiSequence *alignment = ProcessTree (tree, sequences, sparseMatrices, model);
-
- if (enableAlignOrder){
- alignment->SaveOrdering();
- enableAlignOrder = false;
- }
-
- // tree-based refinement
- // TreeBasedBiPartitioning (sparseMatrices, model, alignment, tree);
-
- // iterative refinement
- for (int i = 0; i < numIterativeRefinementReps; i++)
- DoIterativeRefinement (sparseMatrices, model, alignment);
-
- cerr << endl;
-
- // return final alignment
- return alignment;
-}
-
-/////////////////////////////////////////////////////////////////
-// AlignAlignments()
-//
-// Returns the alignment of two MultiSequence objects.
-/////////////////////////////////////////////////////////////////
-
-MultiSequence *AlignAlignments (MultiSequence *align1, MultiSequence *align2,
- const SafeVector<SafeVector<SparseMatrix *> > &sparseMatrices,
- const ProbabilisticModel &model){
-
- // print some info about the alignment
- if (enableVerbose){
- for (int i = 0; i < align1->GetNumSequences(); i++)
- cerr << ((i==0) ? "[" : ",") << align1->GetSequence(i)->GetLabel();
- cerr << "] vs. ";
- for (int i = 0; i < align2->GetNumSequences(); i++)
- cerr << ((i==0) ? "[" : ",") << align2->GetSequence(i)->GetLabel();
- cerr << "]: ";
- }
-
- VF *posterior = model.BuildPosterior (align1, align2, sparseMatrices, cutoff);
- pair<SafeVector<char> *, float> alignment;
-
- // choose the alignment routine depending on the "cosmetic" gap penalties used
- if (gapOpenPenalty == 0 && gapContinuePenalty == 0)
- alignment = model.ComputeAlignment (align1->GetSequence(0)->GetLength(), align2->GetSequence(0)->GetLength(), *posterior);
- else
- alignment = model.ComputeAlignmentWithGapPenalties (align1, align2,
- *posterior, align1->GetNumSequences(), align2->GetNumSequences(),
- gapOpenPenalty, gapContinuePenalty);
-
- delete posterior;
-
- if (enableVerbose){
-
- // compute total length of sequences
- int totLength = 0;
- for (int i = 0; i < align1->GetNumSequences(); i++)
- for (int j = 0; j < align2->GetNumSequences(); j++)
- totLength += min (align1->GetSequence(i)->GetLength(), align2->GetSequence(j)->GetLength());
-
- // give an "accuracy" measure for the alignment
- cerr << alignment.second / totLength << endl;
- }
-
- // now build final alignment
- MultiSequence *result = new MultiSequence();
- for (int i = 0; i < align1->GetNumSequences(); i++)
- result->AddSequence (align1->GetSequence(i)->AddGaps(alignment.first, 'X'));
- for (int i = 0; i < align2->GetNumSequences(); i++)
- result->AddSequence (align2->GetSequence(i)->AddGaps(alignment.first, 'Y'));
- if (!enableAlignOrder)
- result->SortByLabel();
-
- // free temporary alignment
- delete alignment.first;
-
- return result;
-}
-
-/////////////////////////////////////////////////////////////////
-// DoRelaxation()
-//
-// Performs one round of the consistency transformation. The
-// formula used is:
-// 1
-// P'(x[i]-y[j]) = --- sum sum P(x[i]-z[k]) P(z[k]-y[j])
-// |S| z in S k
-//
-// where S = {x, y, all other sequences...}
-//
-/////////////////////////////////////////////////////////////////
-
-SafeVector<SafeVector<SparseMatrix *> > DoRelaxation (MultiSequence *sequences,
- SafeVector<SafeVector<SparseMatrix *> > &sparseMatrices){
- const int numSeqs = sequences->GetNumSequences();
-
- SafeVector<SafeVector<SparseMatrix *> > newSparseMatrices (numSeqs, SafeVector<SparseMatrix *>(numSeqs, NULL));
-
- // for every pair of sequences
- for (int i = 0; i < numSeqs; i++){
- for (int j = i+1; j < numSeqs; j++){
- Sequence *seq1 = sequences->GetSequence (i);
- Sequence *seq2 = sequences->GetSequence (j);
-
- if (enableVerbose)
- cerr << "Relaxing (" << i+1 << ") " << seq1->GetHeader() << " vs. "
- << "(" << j+1 << ") " << seq2->GetHeader() << ": ";
-
- // get the original posterior matrix
- VF *posteriorPtr = sparseMatrices[i][j]->GetPosterior(); assert (posteriorPtr);
- VF &posterior = *posteriorPtr;
-
- const int seq1Length = seq1->GetLength();
- const int seq2Length = seq2->GetLength();
-
- // contribution from the summation where z = x and z = y
- for (int k = 0; k < (seq1Length+1) * (seq2Length+1); k++) posterior[k] += posterior[k];
-
- if (enableVerbose)
- cerr << sparseMatrices[i][j]->GetNumCells() << " --> ";
-
- // contribution from all other sequences
- for (int k = 0; k < numSeqs; k++) if (k != i && k != j){
- if (k < i)
- Relax1 (sparseMatrices[k][i], sparseMatrices[k][j], posterior);
- else if (k > i && k < j)
- Relax (sparseMatrices[i][k], sparseMatrices[k][j], posterior);
- else {
- SparseMatrix *temp = sparseMatrices[j][k]->ComputeTranspose();
- Relax (sparseMatrices[i][k], temp, posterior);
- delete temp;
- }
- }
-
- // now renormalization
- for (int k = 0; k < (seq1Length+1) * (seq2Length+1); k++) posterior[k] /= numSeqs;
-
- // mask out positions not originally in the posterior matrix
- SparseMatrix *matXY = sparseMatrices[i][j];
- for (int y = 0; y <= seq2Length; y++) posterior[y] = 0;
- for (int x = 1; x <= seq1Length; x++){
- SafeVector<PIF>::iterator XYptr = matXY->GetRowPtr(x);
- SafeVector<PIF>::iterator XYend = XYptr + matXY->GetRowSize(x);
- VF::iterator base = posterior.begin() + x * (seq2Length + 1);
- int curr = 0;
- while (XYptr != XYend){
-
- // zero out all cells until the first filled column
- while (curr < XYptr->first){
- base[curr] = 0;
- curr++;
- }
-
- // now, skip over this column
- curr++;
- ++XYptr;
- }
-
- // zero out cells after last column
- while (curr <= seq2Length){
- base[curr] = 0;
- curr++;
- }
- }
-
- // save the new posterior matrix
- newSparseMatrices[i][j] = new SparseMatrix (seq1->GetLength(), seq2->GetLength(), posterior);
- newSparseMatrices[j][i] = NULL;
-
- if (enableVerbose)
- cerr << newSparseMatrices[i][j]->GetNumCells() << " -- ";
-
- delete posteriorPtr;
-
- if (enableVerbose)
- cerr << "done." << endl;
- }
- }
-
- return newSparseMatrices;
-}
-
-/////////////////////////////////////////////////////////////////
-// Relax()
-//
-// Computes the consistency transformation for a single sequence
-// z, and adds the transformed matrix to "posterior".
-/////////////////////////////////////////////////////////////////
-
-void Relax (SparseMatrix *matXZ, SparseMatrix *matZY, VF &posterior){
-
- assert (matXZ);
- assert (matZY);
-
- int lengthX = matXZ->GetSeq1Length();
- int lengthY = matZY->GetSeq2Length();
- assert (matXZ->GetSeq2Length() == matZY->GetSeq1Length());
-
- // for every x[i]
- for (int i = 1; i <= lengthX; i++){
- SafeVector<PIF>::iterator XZptr = matXZ->GetRowPtr(i);
- SafeVector<PIF>::iterator XZend = XZptr + matXZ->GetRowSize(i);
-
- VF::iterator base = posterior.begin() + i * (lengthY + 1);
-
- // iterate through all x[i]-z[k]
- while (XZptr != XZend){
- SafeVector<PIF>::iterator ZYptr = matZY->GetRowPtr(XZptr->first);
- SafeVector<PIF>::iterator ZYend = ZYptr + matZY->GetRowSize(XZptr->first);
- const float XZval = XZptr->second;
-
- // iterate through all z[k]-y[j]
- while (ZYptr != ZYend){
- base[ZYptr->first] += XZval * ZYptr->second;
- ZYptr++;
- }
- XZptr++;
- }
- }
-}
-
-/////////////////////////////////////////////////////////////////
-// Relax1()
-//
-// Computes the consistency transformation for a single sequence
-// z, and adds the transformed matrix to "posterior".
-/////////////////////////////////////////////////////////////////
-
-void Relax1 (SparseMatrix *matZX, SparseMatrix *matZY, VF &posterior){
-
- assert (matZX);
- assert (matZY);
-
- int lengthZ = matZX->GetSeq1Length();
- int lengthY = matZY->GetSeq2Length();
-
- // for every z[k]
- for (int k = 1; k <= lengthZ; k++){
- SafeVector<PIF>::iterator ZXptr = matZX->GetRowPtr(k);
- SafeVector<PIF>::iterator ZXend = ZXptr + matZX->GetRowSize(k);
-
- // iterate through all z[k]-x[i]
- while (ZXptr != ZXend){
- SafeVector<PIF>::iterator ZYptr = matZY->GetRowPtr(k);
- SafeVector<PIF>::iterator ZYend = ZYptr + matZY->GetRowSize(k);
- const float ZXval = ZXptr->second;
- VF::iterator base = posterior.begin() + ZXptr->first * (lengthY + 1);
-
- // iterate through all z[k]-y[j]
- while (ZYptr != ZYend){
- base[ZYptr->first] += ZXval * ZYptr->second;
- ZYptr++;
- }
- ZXptr++;
- }
- }
-}
-
-/////////////////////////////////////////////////////////////////
-// GetSubtree
-//
-// Returns set containing all leaf labels of the current subtree.
-/////////////////////////////////////////////////////////////////
-
-set<int> GetSubtree (const TreeNode *tree){
- set<int> s;
-
- if (tree->GetSequenceLabel() == -1){
- s = GetSubtree (tree->GetLeftChild());
- set<int> t = GetSubtree (tree->GetRightChild());
-
- for (set<int>::iterator iter = t.begin(); iter != t.end(); ++iter)
- s.insert (*iter);
- }
- else {
- s.insert (tree->GetSequenceLabel());
- }
-
- return s;
-}
-
-/////////////////////////////////////////////////////////////////
-// TreeBasedBiPartitioning
-//
-// Uses the iterative refinement scheme from MUSCLE.
-/////////////////////////////////////////////////////////////////
-
-void TreeBasedBiPartitioning (const SafeVector<SafeVector<SparseMatrix *> > &sparseMatrices,
- const ProbabilisticModel &model, MultiSequence* &alignment,
- const TreeNode *tree){
- // check if this is a node of the alignment tree
- if (tree->GetSequenceLabel() == -1){
- TreeBasedBiPartitioning (sparseMatrices, model, alignment, tree->GetLeftChild());
- TreeBasedBiPartitioning (sparseMatrices, model, alignment, tree->GetRightChild());
-
- set<int> leftSubtree = GetSubtree (tree->GetLeftChild());
- set<int> rightSubtree = GetSubtree (tree->GetRightChild());
- set<int> leftSubtreeComplement, rightSubtreeComplement;
-
- // calculate complement of each subtree
- for (int i = 0; i < alignment->GetNumSequences(); i++){
- if (leftSubtree.find(i) == leftSubtree.end()) leftSubtreeComplement.insert (i);
- if (rightSubtree.find(i) == rightSubtree.end()) rightSubtreeComplement.insert (i);
- }
-
- // perform realignments for edge to left child
- if (!leftSubtree.empty() && !leftSubtreeComplement.empty()){
- MultiSequence *groupOneSeqs = alignment->Project (leftSubtree); assert (groupOneSeqs);
- MultiSequence *groupTwoSeqs = alignment->Project (leftSubtreeComplement); assert (groupTwoSeqs);
- delete alignment;
- alignment = AlignAlignments (groupOneSeqs, groupTwoSeqs, sparseMatrices, model);
- }
-
- // perform realignments for edge to right child
- if (!rightSubtree.empty() && !rightSubtreeComplement.empty()){
- MultiSequence *groupOneSeqs = alignment->Project (rightSubtree); assert (groupOneSeqs);
- MultiSequence *groupTwoSeqs = alignment->Project (rightSubtreeComplement); assert (groupTwoSeqs);
- delete alignment;
- alignment = AlignAlignments (groupOneSeqs, groupTwoSeqs, sparseMatrices, model);
- }
- }
-}
-
-/////////////////////////////////////////////////////////////////
-// DoIterativeRefinement()
-//
-// Performs a single round of randomized partionining iterative
-// refinement.
-/////////////////////////////////////////////////////////////////
-
-void DoIterativeRefinement (const SafeVector<SafeVector<SparseMatrix *> > &sparseMatrices,
- const ProbabilisticModel &model, MultiSequence* &alignment){
- set<int> groupOne, groupTwo;
-
- // create two separate groups
- for (int i = 0; i < alignment->GetNumSequences(); i++){
- if (rand() % 2)
- groupOne.insert (i);
- else
- groupTwo.insert (i);
- }
-
- if (groupOne.empty() || groupTwo.empty()) return;
-
- // project into the two groups
- MultiSequence *groupOneSeqs = alignment->Project (groupOne); assert (groupOneSeqs);
- MultiSequence *groupTwoSeqs = alignment->Project (groupTwo); assert (groupTwoSeqs);
- delete alignment;
-
- // realign
- alignment = AlignAlignments (groupOneSeqs, groupTwoSeqs, sparseMatrices, model);
-
- delete groupOneSeqs;
- delete groupTwoSeqs;
-}
-
-/////////////////////////////////////////////////////////////////
-// WriteAnnotation()
-//
-// Computes annotation for multiple alignment and write values
-// to a file.
-/////////////////////////////////////////////////////////////////
-
-void WriteAnnotation (MultiSequence *alignment,
- const SafeVector<SafeVector<SparseMatrix *> > &sparseMatrices){
- ofstream outfile (annotationFilename.c_str());
-
- if (outfile.fail()){
- cerr << "ERROR: Unable to write annotation file." << endl;
- exit (1);
- }
-
- const int alignLength = alignment->GetSequence(0)->GetLength();
- const int numSeqs = alignment->GetNumSequences();
-
- SafeVector<int> position (numSeqs, 0);
- SafeVector<SafeVector<char>::iterator> seqs (numSeqs);
- for (int i = 0; i < numSeqs; i++) seqs[i] = alignment->GetSequence(i)->GetDataPtr();
- SafeVector<pair<int,int> > active;
- active.reserve (numSeqs);
-
- // for every column
- for (int i = 1; i <= alignLength; i++){
-
- // find all aligned residues in this particular column
- active.clear();
- for (int j = 0; j < numSeqs; j++){
- if (seqs[j][i] != '-'){
- active.push_back (make_pair(j, ++position[j]));
- }
- }
-
- outfile << setw(4) << ComputeScore (active, sparseMatrices) << endl;
- }
-
- outfile.close();
-}
-
-/////////////////////////////////////////////////////////////////
-// ComputeScore()
-//
-// Computes the annotation score for a particular column.
-/////////////////////////////////////////////////////////////////
-
-int ComputeScore (const SafeVector<pair<int, int> > &active,
- const SafeVector<SafeVector<SparseMatrix *> > &sparseMatrices){
-
- if (active.size() <= 1) return 0;
-
- // ALTERNATIVE #1: Compute the average alignment score.
-
- float val = 0;
- for (int i = 0; i < (int) active.size(); i++){
- for (int j = i+1; j < (int) active.size(); j++){
- val += sparseMatrices[active[i].first][active[j].first]->GetValue(active[i].second, active[j].second);
- }
- }
-
- return (int) (200 * val / ((int) active.size() * ((int) active.size() - 1)));
-
-}
Modified: trunk/packages/probcons/trunk/debian/README.Debian
===================================================================
--- trunk/packages/probcons/trunk/debian/README.Debian 2007-04-12 10:59:18 UTC (rev 269)
+++ trunk/packages/probcons/trunk/debian/README.Debian 2007-04-18 13:39:18 UTC (rev 270)
@@ -9,6 +9,7 @@
distribution of ProbCons: probcons-RNA. It is an experimental version of
ProbCons compiled with parameters estimated via unsupervised training on
BRAliBASE. The parameters were downloaded from a separate archive:
- http://probcons.stanford.edu/probconsRNA.tar.gz
+ http://probcons.stanford.edu/probconsRNA.tar.gz, and the current version of
+ probcons was used to build probcons-RNA.
- -- Charles Plessy <charles-debian-nospam at plessy.org>, Wed, 3 May 2006 09:20:57 +0900
+ -- Charles Plessy <charles-debian-nospam at plessy.org>, Wed, 18 Apr 2007 22:30:52 +0900
Modified: trunk/packages/probcons/trunk/debian/changelog
===================================================================
--- trunk/packages/probcons/trunk/debian/changelog 2007-04-12 10:59:18 UTC (rev 269)
+++ trunk/packages/probcons/trunk/debian/changelog 2007-04-18 13:39:18 UTC (rev 270)
@@ -1,3 +1,11 @@
+probcons (1.12-2) unstable; urgency=low
+
+ * Corrected an error in the build rules making probcons-RNA, which
+ prevented it to benefit from the fixes on probcons.
+ (Closes: #417482) (again!)
+
+ -- Charles Plessy <charles-debian-nospam at plessy.org> Wed, 18 Apr 2007 22:34:51 +0900
+
probcons (1.12-1) unstable; urgency=low
* New upstream release, buildable with gcc 4.3 (Closes: #417482)
Modified: trunk/packages/probcons/trunk/debian/rules
===================================================================
--- trunk/packages/probcons/trunk/debian/rules 2007-04-12 10:59:18 UTC (rev 269)
+++ trunk/packages/probcons/trunk/debian/rules 2007-04-18 13:39:18 UTC (rev 270)
@@ -28,7 +28,7 @@
build: patch manpages
dh_testdir
- cp debian/Main-RNA.cc Main-RNA.cc
+ cp Main.cc Main-RNA.cc
cp debian/Defaults-RNA.h Defaults-RNA.h
$(MAKE)
touch build-stamp
More information about the debian-med-commit
mailing list