main.cpp

ViennaRNA-1.6.1

	      {
		it->second->printFastaAli(false);
		cout << endl;
	      }
	    else
	      {  

		it->second->printSeqAli();

#ifdef HAVE_LIBRNA  // This features require the ViennaRNA library		
		// print alignment
		it->second->printStrAli();
		cout << endl;
#endif		
	      }

	    // save profile
	    if(options.has(RNAforesterOptions::SaveProfile))
	    {
	      string filename;
	      filename=options.generateFilename(RNAforesterOptions::SaveProfile,".pro", "rna.pro",clusterNr);
	      it->second->save(filename);
	    }
	    

		// print consensus structure
		cout << "Consensus sequence/structure:" << endl;
		it->second->printConsensus(minPairProb);
		cout << endl;

#ifdef HAVE_LIBG2  // This features require the g2 library		
		// generate squiggle plots
			if(options.has(RNAforesterOptions::MakeSquigglePlot))
			{
				RNAProfileAlignment::SquigglePlotOptions sqOptions;
				string filename;

				// plot showing full base information
				filename=options.generateFilename(RNAforesterOptions::MakeSquigglePlot,".ps", "rnaprofile.ps",clusterNr);
				sqOptions.greyColors=options.has(RNAforesterOptions::SquiggleGreyColors);
				sqOptions.minPairProb=minPairProb;
				sqOptions.mostLikelySequence=false;
				it->second->squigglePlot(filename,sqOptions);

				// plot showing consensus sequence
				filename=options.generateFilename(RNAforesterOptions::MakeSquigglePlot,"_cons.ps", "rnaprofile_cs.ps",clusterNr);
				sqOptions.mostLikelySequence=true;
				it->second->squigglePlot(filename,sqOptions);
			}
#endif		      

		clusterNr++;
	  }


	/*
	if(!options.has(RNAforesterOptions::ShowOnlyScore))
	{
		cout << "global optimal score: ";
	}
	cout << optScore << endl;

	list<RNAProfileAlignment*>::iterator it=inputMapProfile.begin();
	RNAProfileAlignment *ppf=it;
	if(!options.has(RNAforesterOptions::ShowOnlyScore))
	{
		// generate dot file
		makeDotFileAli(*ppf,options);

		ppf->print();
		cout << endl;
		ppf->printConsensus();
	}
	*/

	// save profile alignment to binary file
	/*	      if(options.has(RNAforesterOptions::SaveMultipleAliFile))
	{
	string filename;

	filename=generateFilename(options,RNAforesterOptions::SaveMultipleAliFile,".mta", "unknown.dot");
	ofstream s(filename.c_str());
	f1->save(s);
	}
	*/

	/*
	// generate squiggle plots
	if(options.has(RNAforesterOptions::MakeSquigglePlot))
	{
		RNAProfileAlignment::SquigglePlotOptions sqOptions;
		string filename;

		// plot showing full base information
		filename=options.generateFilename(RNAforesterOptions::MakeSquigglePlot,".ps", "rnaprofile.ps");
		sqOptions.greyColors=options.has(RNAforesterOptions::SquiggleGreyColors);
		sqOptions.mostLikelySequence=false;
		ppf->squigglePlot(filename,sqOptions);

		// plot showing consensus sequence
		filename=options.generateFilename(RNAforesterOptions::MakeSquigglePlot,"_cons.ps", "rnaprofile_cons.ps");
		sqOptions.greyColors=options.has(RNAforesterOptions::SquiggleGreyColors);
		sqOptions.mostLikelySequence=true;
		ppf->squigglePlot(filename,sqOptions);
	}
	*/

	
	delete alg;
}


void alignPairwise(deque<RNAForest*> &inputListPW,Score &score,const RNAforesterOptions &options)
{
	IntScoreRNA_AlgebraType *alg;
	Uint xbasepos,ybasepos,xlen,ylen;
	list<pair<Uint,Uint> > xsubopts;
	list<pair<Uint,Uint> > ysubopts;
	string seq1,seq2,str1,str2;
	char s[8];
	int suboptPercent;
	Uint count=1;
	RNAFuncs::SquigglePlotOptions sqOptions;
	tms tmsStart, tmsEnd;

	// read options
	options.get(RNAforesterOptions::LocalSubopts,suboptPercent,100);

#ifdef HAVE_LIBG2  // This features require the g2 library	
	// generate squiggle plot
	if(options.has(RNAforesterOptions::MakeSquigglePlot))
	{	      
		// get sq options
		sqOptions.hideBaseNumbers=options.has(RNAforesterOptions::SquiggleHideBaseNumbers);
		options.get(RNAforesterOptions::SquiggleBaseNumberInterval,sqOptions.baseNumInterval,(Uint)20);
		sqOptions.greyColors=options.has(RNAforesterOptions::SquiggleGreyColors);
		options.get(RNAforesterOptions::SquiggleScaleFactor,sqOptions.scale,1.0);
		sqOptions.generateFIG=options.has(RNAforesterOptions::SquiggleGenerateFIG);		
#ifdef HAVE_LIBGD
		sqOptions.generatePNG=options.has(RNAforesterOptions::SquiggleGeneratePNG);
		sqOptions.generateJPG=options.has(RNAforesterOptions::SquiggleGenerateJPG);
#endif
		
	}
#endif

	// distance or similarity
	if(options.has(RNAforesterOptions::CalculateDistance))
		alg=new IntDistRNA_Algebra(score);
	else
	{
		if(options.has(RNAforesterOptions::RIBOSUMScore))
			alg=new RIBOSUM8560(score);
		else
			alg=new IntSimiRNA_Algebra(score);
	}

	RNAForest *f1=inputListPW.front();
	RNAForest *f2=inputListPW.back();

	if(options.has(RNAforesterOptions::SpaceTimeInfo))
	  {
	    IntScore_AlgebraType *algGlobClassic;

	    algGlobClassic=new IntDist_Algebra(score);

	    if(!options.has(RNAforesterOptions::ShowOnlyScore))
	      {
		cout << "F1_NUMNODES" << ";";
		cout << "F2_NUMNODES" << ";";
		cout << "F1_DEGREE" << ";";
		cout << "F2_DEGREE" << ";";
		cout << "F1_LEAVES" << ";"; 
		cout << "F2_LEAVES" << ";";
		cout << "F1_DEPTH" << ";"; 
		cout << "F2_DEPTH" << ";";
		cout << "F1_NUMCSFS" << ";";
		cout << "F2_NUMCSFS" << ";";
		cout << "TABLE_SIZE_4D" << ";";
		cout << "TABLE_SIZE_2D" << ";";
		cout << "GLOBALI_TIME" << ";";
		cout << "GLOBALI_TIME_SPEEDUP" << ";";
		cout << "LOCALALI_TIME" << ";";
		cout << "LOCALALI_TIME_SPEEDUP" << ";";
		
		cout << endl;
	      }

	    cout << f1->size() << "\t";
	    cout << f2->size() << "\t"; 
	    cout << f1->maxDegree() << "\t";
	    cout << f2->maxDegree() << "\t";
	    cout << f1->numLeaves() << "\t";
            cout << f2->numLeaves() << "\t";
	    cout << f1->maxDepth() << "\t";
	    cout << f2->maxDepth() << "\t";
	    cout << f1->getNumCSFs() << "\t";
	    cout << f2->getNumCSFs() << "\t";
	    cout << f1->size()*f2->size()*f1->maxDegree()*f1->maxDegree() << "\t";
	    cout << f1->getNumCSFs()*f2->getNumCSFs() << "\t";

 	    // global alignment
	    {
	      times(&tmsStart);
	      Alignment<int,RNA_Alphabet,RNA_AlphaPair> ali(f1,f2,*algGlobClassic,false,true);
	      times(&tmsEnd);
	      cout <<((double) (tmsEnd.tms_utime - tmsStart.tms_utime))/sysconf(_SC_CLK_TCK) << "\t"; 
	      //	    cerr << "#" << ali.getGlobalOptimum() << "#";
	    }

 	    // global alignment speedup
	    {
	      times(&tmsStart);
	      Alignment<int,RNA_Alphabet,RNA_AlphaPair> ali(f1,f2,*algGlobClassic,false,false);
	      times(&tmsEnd);
	      cout <<((double) (tmsEnd.tms_utime - tmsStart.tms_utime))/sysconf(_SC_CLK_TCK) << "\t"; 
	      //	    cerr << "#" << ali.getGlobalOptimum() << "#";
	    }

 	    // local alignment
	    {
	      times(&tmsStart);
	      Alignment<int,RNA_Alphabet,RNA_AlphaPair> ali(f1,f2,*algGlobClassic,true,true);
	      times(&tmsEnd);
	      cout <<((double) (tmsEnd.tms_utime - tmsStart.tms_utime))/sysconf(_SC_CLK_TCK) << "\t"; 
	      //	    cerr << "#" << ali.getLocalOptimum() << "#";
	    }

 	    // local alignment speedup
	    {
	      times(&tmsStart);
	      Alignment<int,RNA_Alphabet,RNA_AlphaPair> ali(f1,f2,*algGlobClassic,true,false);
	      times(&tmsEnd);
	      cout <<((double) (tmsEnd.tms_utime - tmsStart.tms_utime))/sysconf(_SC_CLK_TCK) << "\t"; 
	      //	    cerr << "#" << ali.getLocalOptimum() << "#";
	    }
	    

	    cout << endl;

	    exit(EXIT_SUCCESS);
	  }

	Alignment<int,RNA_Alphabet,RNA_AlphaPair> ali(f1,f2,*alg);
	RNA_Alignment ppfali;
	ppfali.setStructureNames(f1->getName(),f2->getName());


	if(!options.has(RNAforesterOptions::ShowOnlyScore))
	{
	  if(options.has(RNAforesterOptions::SmallInLarge))
	    {
	      cout << "small-in-large optimal score: ";
	    }
	  else
	    {
	      if(options.has(RNAforesterOptions::LocalSimilarity))
		cout << "local optimal score: ";
	      else
		cout << "global optimal score: ";
	    }
	}

	if(options.has(RNAforesterOptions::SmallInLarge))
	  {
	    cout << ali.getSILOptimum() << endl;
	  }
	else
	  {
	    if(options.has(RNAforesterOptions::LocalSimilarity))
	      {
		cout << ali.getLocalOptimum() << endl;
	      }  
	    else
	      {
		cout << ali.getGlobalOptimum() << endl;	  	  
		if(options.has(RNAforesterOptions::RelativeScore))
		  cout << ali.getGlobalOptimumRelative() << endl;
	      }
	  }


	if(!options.has(RNAforesterOptions::ShowOnlyScore))
	{
	  if(options.has(RNAforesterOptions::SmallInLarge))
	    {
	      ali.getOptSILAlignment(ppfali,ybasepos);
	      cout << "starting at position: " << ybasepos << endl << endl; 
	    }
	  else
	    {	      
	      if(options.has(RNAforesterOptions::LocalSimilarity))
		{
		  ali.resetOptLocalAlignment(suboptPercent);
		  ali.getOptLocalAlignment(ppfali,xbasepos,ybasepos);		      		    
		  cout << "starting at positions: " << xbasepos << "," << ybasepos << endl << endl; 
		}
	      else		  
		ali.getOptGlobalAlignment(ppfali);		  
	    }

		// generate dot file
		makeDotFileAli(ppfali,options);

		// print alignment
		ppfali.getSequenceAlignments(seq1,seq2);				
		ppfali.getStructureAlignment(str1,true);
		ppfali.getStructureAlignment(str2,false);

		if(options.has(RNAforesterOptions::FastaOutput))
		  {
		    cout << ppfali.getStructureNameX() << endl;
		    cout << seq1 << endl;
		    cout << str1 << endl;
		    cout << ppfali.getStructureNameY() << endl;
		    cout << seq2 << endl;
		    cout << str2 << endl;
		    cout << endl;
		  }
		else		
		  RNAFuncs::printAli(ppfali.getStructureNameX(),ppfali.getStructureNameY(),seq1,seq2,str1,str2);  


		xlen=seq1.size();
		ylen=seq2.size();
		if(options.has(RNAforesterOptions::LocalSimilarity))
		{
			xsubopts.push_back(pair<Uint,Uint>(xbasepos,xlen));
			ysubopts.push_back(pair<Uint,Uint>(ybasepos,ylen));
		}

#ifdef HAVE_LIBG2  // This features require the g2 library		
		if(options.has(RNAforesterOptions::MakeSquigglePlot))
		  {
		    sprintf(s,"%d",count);
		    ppfali.squigglePlot(s,sqOptions);      
		  }
#endif	       

		// suboptimal alignments
		if(options.has(RNAforesterOptions::LocalSubopts))
		{
			while(ali.nextLocalSuboptimum())
			{
				count++;
				cout << "local optimal score: ";
				cout << ali.getLocalOptimum() << endl;
				ali.getOptLocalAlignment(ppfali,xbasepos,ybasepos);
				cout << "starting at positions: " << xbasepos << "," << ybasepos << endl << endl; 

				// print alignment
				ppfali.getSequenceAlignments(seq1,seq2);
				ppfali.getStructureAlignment(str1,true);
				ppfali.getStructureAlignment(str2,false);
				RNAFuncs::printAli(ppfali.getStructureNameX(),ppfali.getStructureNameY(),seq1,seq2,str1,str2);  
				xlen=seq1.size();
				ylen=seq2.size();
				xsubopts.push_back(pair<Uint,Uint>(xbasepos,xlen));
				ysubopts.push_back(pair<Uint,Uint>(ybasepos,ylen));	

#ifdef HAVE_LIBG2  // This features require the g2 library		
		                if(options.has(RNAforesterOptions::MakeSquigglePlot))		  
				{
				  sprintf(s,"%d",count);				
				  ppfali.squigglePlot(s,sqOptions);      	    	  
				}
#endif			       
			}
		}

	}

#ifdef HAVE_LIBRNA  // This features require the ViennaRNA library	
	// generate xml
	if(options.has(RNAforesterOptions::GenerateXML))
	{
		string filename;

		filename=options.generateFilename(RNAforesterOptions::GenerateXML,".xml", "ali.xml");
		ofstream s(filename.c_str());
		ppfali.generateXML(s);
	}
#endif	

#ifdef HAVE_LIBG2  // This features require the g2 library
	// generate squiggle plot
	if(options.has(RNAforesterOptions::MakeSquigglePlot))
	{
		f1->plot2d("x_str", xsubopts, sqOptions);
		f2->plot2d("y_str", ysubopts, sqOptions);  
	}
#endif	

	// clear input list
	deque<RNAForest*>::const_iterator it;
	for(it = inputListPW.begin(); it!=inputListPW.end(); it++)
		delete *it;
	inputListPW.clear();
	delete alg;
}


void editPairwise(list<RNAForestSZ*> &inputListSZ,Score &score,RNAforesterOptions &options)
{
  //  timeb t1,t2;
  IntScoreSZAlgebraType *alg;

//  if(options.has(RNAforesterOptions::CalculateDistance))
    alg=new IntDistSZAlgebra(score);
//  else
//    alg=new IntSimiSZAlgebra(score);
  
  RNAForestSZ *f1=inputListSZ.front();
  RNAForestSZ *f2=inputListSZ.back();

  //  ftime(&t1);
  Mapping<int,RNA_Alphabet> mapping(f1,f2,*alg);
  //  ftime(&t2);

//	f1->printParameters("F1");
//	f2->printParameters("F2");

  cout << "Global optimum: " << mapping.getGlobalOptimum() << endl;
  //  cout << "Calculation Time ms: " << (t2.time*1000+t2.millitm) - (t1.time*1000+t1.millitm) << endl;
}



void alignPairwiseSimple(deque<RNAForest*> &inputListPW,Score &score,RNAforesterOptions &options)
{
  //  timeb t1,t2;
  IntScore_AlgebraType *alg;

//  if(options.has(RNAforesterOptions::CalculateDistance))
    alg=new IntDist_Algebra(score);
//  else
//    alg=new ScoreAlgebraSimple(score);

  RNAForest *f1=inputListPW.front();
  RNAForest *f2=inputListPW.back();

  //  ftime(&t1);
  Alignment<int,RNA_Alphabet,RNA_AlphaPair> ali(f1,f2,*alg,false);
  //  ftime(&t2);

//	f1->printParameters("F1");
  //f2->printParameters("F2");

  cout << "Global optimum: " << ali.getGlobalOptimum() << endl;
  //  cout << "Calculation Time ms: " << (t2.time*1000+t2.millitm) - (t1.time*1000+t1.millitm) << endl;
}