--- /dev/null
+/*
+ Copyright by Matthias Hoechsmann (C) 2002-2004
+ =====================================
+ You may use, copy and distribute this file freely as long as you
+ - do not change the file,
+ - leave this copyright notice in the file,
+ - do not make any profit with the distribution of this file
+ - give credit where credit is due
+ You are not allowed to copy or distribute this file otherwise
+ The commercial usage and distribution of this file is prohibited
+ Please report bugs and suggestions to <mhoechsm@TechFak.Uni-Bielefeld.DE>
+*/
+
+#include <algorithm>
+#include <fstream>
+
+#include "alignment.h"
+#include "progressive_align.h"
+
+#include "alignment.t.cpp"
+
+using namespace std;
+
+// !!! this operator is defined as > !!!
+bool operator < (pair<double,RNAProfileAlignment*> &l, pair<double,RNAProfileAlignment*> &r)
+{
+ if(l.second->getNumStructures() < r.second->getNumStructures())
+ return false;
+ else
+ return true;
+}
+
+void progressiveAlign(deque<RNAProfileAlignment*> &inputList, deque<pair<double,RNAProfileAlignment*> > &resultList, const DoubleScoreProfileAlgebraType *alg,const RNAforesterOptions &options)
+{
+ RNAProfileAliMapType inputMapProfile;
+ deque<RNAProfileAliKeyPairType> alignList;
+ RNAProfileAliMapType::iterator it,it2;
+ long x,y,bestx,besty; // keys of stuctures in inputMapProfile
+ RNAProfileAlignment *ppfx=NULL,*ppfy=NULL,*ppf=NULL;
+ int level=1;
+ double bestScore,threshold,cutoff;
+ string clusterfilename;
+ ofstream s;
+ Matrix<double> *score_mtrx;
+ bool local;
+
+ local=options.has(RNAforesterOptions::LocalSimilarity);
+
+ cout << "*** Calculation ***" << endl << endl;
+
+ // create inputMapProfile to access a profile by an index value
+ deque<RNAProfileAlignment*>::const_iterator inpIt;
+ long i=1;
+ for(inpIt=inputList.begin();inpIt!=inputList.end();inpIt++)
+ {
+ inputMapProfile[i]=*inpIt;
+ i++;
+ }
+ inputList.clear();
+
+ // create matrix for all against all comparison
+ bestScore=alg->worst_score();
+ score_mtrx=new Matrix<double>(inputMapProfile.size(),inputMapProfile.size());
+
+ // set threshold for the clustering algorithm
+ if(options.has(RNAforesterOptions::CalculateDistance))
+ options.get(RNAforesterOptions::ClusterThreshold,threshold,20.0);
+ else
+ options.get(RNAforesterOptions::ClusterThreshold,threshold,0.7);
+ cout << "clustering threshold is: " << threshold << endl;
+
+ // set cutoff value for clustering
+ if(options.has(RNAforesterOptions::CalculateDistance))
+ options.get(RNAforesterOptions::ClusterJoinCutoff,cutoff,100.0);
+ else
+ options.get(RNAforesterOptions::ClusterJoinCutoff,cutoff,0.0);
+ cout << "join clusters cutoff is: " << cutoff << endl << endl;
+
+ // generate dot file
+ clusterfilename=options.generateFilename(RNAforesterOptions::Help,"_cluster.dot", "cluster.dot"); // use Help as dummy
+ s.open(clusterfilename.c_str());
+ s << "digraph forest" << endl << "{" << endl;
+
+ // generate nodes for the input forests
+ for(it=inputMapProfile.begin();it!=inputMapProfile.end();it++)
+ {
+ ppf=it->second;
+
+ s << "\"" << ppf->getName() << "\"" << "[label=\"" << ppf->getName() << "\"]" << endl;
+ s << "\"" << ppf->getName() << "\"" << "[label=\"" << ppf->getName() << "\"]" << endl;
+ }
+
+
+ // compute all pairwise distances
+ // !! NOTE !! iterating through the map is ordered by key number
+ // as i only calculate a triangle matrix this is a prerequisite
+ cout << "Computing all pairwise similarities" << endl;
+ for(it=inputMapProfile.begin();it!=inputMapProfile.end();it++)
+ {
+ x=it->first;
+ ppfx=it->second;
+ for(it2=inputMapProfile.begin();it2->first<it->first;it2++)
+ {
+ y=it2->first;
+ ppfy=it2->second;
+
+ Alignment<double,RNA_Alphabet_Profile,RNA_Alphabet_Profile> ali(ppfx,ppfy,*alg,local);
+ if(local)
+ score_mtrx->setAt(x-1,y-1,ali.getLocalOptimum());
+ else
+ score_mtrx->setAt(x-1,y-1,ali.getGlobalOptimumRelative());
+
+ cout << x << "," << y << ": " << score_mtrx->getAt(x-1,y-1) << endl;
+
+ WATCH(DBG_MULTIPLE,"progressiveAlign",x);
+ WATCH(DBG_MULTIPLE,"progressiveAlign",y);
+ WATCH(DBG_MULTIPLE,"progressiveAlign",ali.getGlobalOptimumRelative());
+ }
+ }
+ cout << endl;
+
+ while(inputMapProfile.size()>1)
+ {
+ // find the best score of all pairwise alignments
+ bestScore=alg->worst_score();
+ for(it=inputMapProfile.begin();it!=inputMapProfile.end();it++)
+ {
+ x=it->first;
+ for(it2=inputMapProfile.begin();it2->first < it->first;it2++)
+ {
+ double old_bestScore=bestScore;
+
+ y=it2->first;
+ WATCH(DBG_MULTIPLE,"progressiveAlign",x);
+ WATCH(DBG_MULTIPLE,"progressiveAlign",y);
+ WATCH(DBG_MULTIPLE,"progressiveAlign",score_mtrx->getAt(x-1,y-1));
+
+ WATCH(DBG_MULTIPLE,"progressiveAlign",bestScore);
+ bestScore=alg->choice(bestScore,score_mtrx->getAt(x-1,y-1));
+
+ if(bestScore!=old_bestScore)
+ {
+ bestx=it->first;
+ besty=it2->first;
+ old_bestScore=bestScore;
+ }
+ }
+ }
+
+ WATCH(DBG_MULTIPLE,"progressiveAlign",bestScore);
+ cout << "joining alignments:" << endl;
+
+ // if threshold is set generate a list of best pairs within threshold
+ if(alg->choice(bestScore,threshold)!= threshold)
+ {
+ Graph graph;
+ int *mate;
+ int maximize;
+
+ graph=makePairsGraph(inputMapProfile,alg,score_mtrx,threshold);
+ if(options.has(RNAforesterOptions::CalculateDistance))
+ maximize=0;
+ else
+ maximize=1;
+
+ mate = Weighted_Match(graph,1,maximize);
+ for(x=1;x<=score_mtrx->xDim();x++) // !! begins at 1 !!
+ {
+ if(x<mate[x])
+ {
+ // if it is a best pair put it in the align list
+ ppfx=inputMapProfile[x];
+ inputMapProfile.erase(x);
+ alignList.push_back(make_pair(x,ppfx));
+ ppfy=inputMapProfile[mate[x]];
+ inputMapProfile.erase(mate[x]);
+ alignList.push_back(make_pair(mate[x],ppfy));
+ }
+ }
+ free(mate);
+ free(graph);
+ }
+ else
+ {
+ // if there us no pair below the threshold
+ // combine those two profile forests, that produced the best score
+ ppfx=inputMapProfile[bestx];
+ ppfy=inputMapProfile[besty];
+ inputMapProfile.erase(bestx);
+ inputMapProfile.erase(besty);
+ alignList.push_back(make_pair(bestx,ppfx));
+ alignList.push_back(make_pair(besty,ppfy));
+ }
+
+ // align all forests in the align list.
+ while(alignList.size()>1)
+ {
+ string aliName;
+
+ x=alignList.front().first;
+ ppfx=alignList.front().second;
+ alignList.erase(alignList.begin());
+ y=alignList.front().first;
+ ppfy=alignList.front().second;
+ alignList.erase(alignList.begin());
+
+ // compute alignment again
+ Alignment<double,RNA_Alphabet_Profile,RNA_Alphabet_Profile> bestali(ppfx,ppfy,*alg,local);
+ if(local)
+ bestScore=bestali.getLocalOptimum();
+ else
+ bestScore=bestali.getGlobalOptimumRelative();
+
+ // test, if score is worse than the cutoff value
+ if(alg->choice(bestScore,cutoff)== cutoff)
+ {
+ // copy involved forests to the result list
+ cout << x << "," << y << ": alignment is below cutoff." << endl;
+ if(ppfx->getNumStructures()>1)
+ resultList.push_back(make_pair(ppfx->maxScore(*alg),ppfx));
+ if(ppfy->getNumStructures()>1)
+ resultList.push_back(make_pair(ppfy->maxScore(*alg),ppfy));
+ }
+ else
+ {
+ // calculate optimal alignment and add it to inputMapProfile
+
+ ppf=new RNAProfileAlignment(ppfx->getNumStructures(),ppfy->getNumStructures());
+
+ if(local)
+ {
+ Uint xbasepos,ybasepos;
+ bestali.getOptLocalAlignment(*ppf,xbasepos,ybasepos);
+ }
+ else
+ {
+ bestali.getOptGlobalAlignment(*ppf);
+ }
+
+ ppf->addStrNames(ppfx->getStrNames());
+ ppf->addStrNames(ppfy->getStrNames());
+ aliName=ppfx->getName() + "." +ppfy->getName();
+ ppf->setName(aliName);
+
+ // for debug purposes
+ //string dotfilename;
+ //dotfilename=ppf->getName() + string(".dot");
+ // ofstream s("test.dot");
+ // ppf->printDot(s);
+
+ // generate nodes in cluster file (dot format)
+ s << "\"" << ppf->getName() << "\"" << "[shape=\"diamond\" label=\"" << bestScore << "\"]" << endl;
+ s << "\"" << ppf->getName() << "\"" << "-> {\"" << ppfx->getName() << "\" \"" << ppfy->getName() << "\"}";
+
+ cout << x << "," << y << ": " << bestScore << " -> " << min(x,y) << endl;
+
+ delete ppfx;
+ delete ppfy;
+
+ // calculate distance to all forests in the list
+ cout << "Calculate similarities to other clusters" << endl;
+ ppfx=ppf;
+ // x remains x !!
+ for(it=inputMapProfile.begin();it!=inputMapProfile.end();it++)
+ {
+ y=it->first;
+ ppfy=it->second;
+ Alignment<double,RNA_Alphabet_Profile,RNA_Alphabet_Profile> ali(ppfx,ppfy,*alg,local);
+ if(local)
+ {
+ score_mtrx->setAt(min(x-1,y-1),max(x-1,y-1),ali.getLocalOptimum()); // min - max = fill the upper triangle
+ cout << min(x,y) << "," << max(x,y) << ": " << ali.getLocalOptimum() << endl;
+ }
+ else
+ {
+ score_mtrx->setAt(min(x-1,y-1),max(x-1,y-1),ali.getGlobalOptimumRelative()); // min - max = fill the upper triangle
+ cout << min(x,y) << "," << max(x,y) << ": " << ali.getGlobalOptimumRelative() << endl;
+ }
+ }
+ cout << endl;
+
+ // ... and append it to the list
+ inputMapProfile.insert(make_pair(x,ppf));
+ }
+ }
+
+ level++;
+ }
+
+ assert(inputMapProfile.size()<2);
+
+ // copy last profile to resultList
+ if(inputMapProfile.size()==1)
+ {
+ ppf=inputMapProfile.begin()->second;
+ resultList.push_back(make_pair(ppf->maxScore(*alg),ppf));
+ inputMapProfile.clear();
+ }
+
+ // end of dot file
+ s << "}" << endl;
+
+ // sort result list
+ // std::sort(resultList.begin(),resultList.end());
+
+
+ delete score_mtrx;
+}
+
+Graph makePairsGraph(const RNAProfileAliMapType &inputMapProfile, const DoubleScoreProfileAlgebraType *alg, const Matrix<double> *score_mtrx, double threshold)
+{
+ Graph graph;
+ RNAProfileAliMapType::const_iterator it,it2;
+ RNAProfileAlignment *ppfx=NULL,*ppfy=NULL;
+
+ graph = NewGraph(score_mtrx->xDim());
+
+ for(int i=0;i<score_mtrx->xDim();i++)
+ {
+ Xcoord(graph,i+1) = 0;
+ Ycoord(graph,i+1) = 0;
+ }
+
+ for(it=inputMapProfile.begin();it!=inputMapProfile.end();it++)
+ {
+ ppfx=it->second;
+ for(it2=inputMapProfile.begin();it2->first<it->first;it2++)
+ {
+ double score;
+ ppfy=it2->second;
+
+ score=score_mtrx->getAt(it->first-1,it2->first-1);
+ if(alg->choice(score,threshold) != threshold) // is it better than the threshold ?
+ {
+ AddEdge (graph,it->first,it2->first,(int)(score*100.0));
+ }
+ }
+ }
+
+
+ WriteGraph (graph,"test.out");
+ return graph;
+}
git://source.jalview.org / jabaws.git / blobdiff