main.cpp

ViennaRNA-1.6.1

/*
  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 <[e-mail]>
*/

#include <deque>
#include <functional>
#include <fstream>
#include <iomanip>
#include <iostream>
#include <list>
#include <sstream>
#include <string>
#include <map>

//#include <sys/timeb.h>
#include <sys/times.h>
#include <unistd.h>

#ifndef WIN32
#include "config.h"
#endif

#include "Arguments.h"
#include "alignment.h"
#include "debug.h"

//#include "global_alignment.h"
#include "treeedit.h"

#include "misc.h"
#include "progressive_align.h"
#include "rna_alignment.h"
#include "rnaforest.h"
#include "rnaforestsz.h"
#include "rnafuncs.h"
#include "rna_profile_alignment.h"
#include "rna_algebra.h"
#include "rnaforester_options.h"


#include "alignment.t.cpp"
//#include "global_alignment.t.cpp"
#include "treeedit.t.cpp"
//#include "ppforest.t.cpp"

using namespace std;

/* ****************************************** */
/*          Definitions and typedefs          */
/* ****************************************** */

struct ToLower : public unary_function<char,char> {
	char operator()(char a)
	{
		return tolower(a);
	};
};

template <class L>
void makeDotFileAli(const PPForest<L> &ppf, const RNAforesterOptions &options)
{
	if(options.has(RNAforesterOptions::OutputAlignmentDotFormat))
	{
		string filename;
		options.get(RNAforesterOptions::OutputAlignmentDotFormat,filename,string("ali.dot"));
		ofstream s(filename.c_str());
		ppf.printDot(s);
	}
}

template <class L>
void makeDotFileInp(const PPForest<L> &ppf, const RNAforesterOptions &options, Uint count)
{
	if(options.has(RNAforesterOptions::MakeDotForInputTrees))
	{
		ostringstream ss;
		ofstream os;
		ss << "input" << count << ".dot";
		os.open(ss.str().c_str());
		ppf.printDot(os);
		os.close();
	}
}

static const string RNAFORESTER_VERSION = "1.5";
static const string PROMPT = "Input string (upper or lower case); & to end for multiple alignments, @ to quit\n";
static const string SCALE = "....,....1....,....2....,....3....,....4....,....5....,....6....,....7....,....8\n";

void alignMultiple(deque<RNAProfileAlignment*> &alignList, Score &score,const RNAforesterOptions &options);
void alignPairwise(deque<RNAForest*> &inputListPW,Score &score,const RNAforesterOptions &options);
void cutAfterChar(string &s,char c);

void editPairwise(list<RNAForestSZ*> &inputListSZ,Score &score,RNAforesterOptions &options);
void alignPairwiseSimple(deque<RNAForest*> &inputListPW,Score &score,RNAforesterOptions &options);

static void showversion(const char *prog)
{
	cout << prog << ", version " << RNAFORESTER_VERSION << endl;
	cout << "Copyright Matthias Hoechsmann 2001-2004," << endl << "[e-mail]" << endl;
}

int main(int argc, const char **argv)
{	
	string buffer;
	string baseStr,viennaStr,nameStr;
	Ulong basePairCount,maxDepth;
	deque<RNAForest*> inputListPW;
	deque<RNAProfileAlignment*> alignList;
	bool showScale=true,multipleAlign=false;
	istream *inputStream=NULL;
	ifstream *inputFile=NULL;
	Uint structure_count=1;
	int suboptPercent=100;
	double minPairProb=0.25;

	list<RNAForestSZ*> inputListSZ;

	try
	{
		RNAforesterOptions options(argc,argv);

		// check options
		if(options.has(RNAforesterOptions::Help))
		{
			options.help();
			exit(EXIT_SUCCESS);
		}

		if(options.has(RNAforesterOptions::SecretHelp))
		{
			options.secretHelp();
			exit(EXIT_SUCCESS);
		}

		if(options.has(RNAforesterOptions::Version))
		{
			showversion(argv[0]);
			exit(EXIT_SUCCESS);      
		}

		// read score values
		Score score(options);
		if(!options.has(RNAforesterOptions::ShowOnlyScore))
			score.print();

		// check option suboptimals
		if(options.has(RNAforesterOptions::LocalSubopts))
		{
			options.get(RNAforesterOptions::LocalSubopts,suboptPercent,100);
			if(suboptPercent<0 || suboptPercent>100)
			{
				cerr << "error: value for parameter --subopts must be in range from 0 to 100" << endl;
				exit(EXIT_FAILURE);
			}
		}

#ifdef HAVE_LIBRNA  // This features require the ViennaRNA library		
		options.get(RNAforesterOptions::PredictMinPairProb,minPairProb,0.25);
		if(options.has(RNAforesterOptions::PredictProfile))
		  cout << "Minumim required basepair probability (-pmin): " << minPairProb << endl;
#endif
		
		// show if suboptimals
		if(!options.has(RNAforesterOptions::ShowOnlyScore))
			if(options.has(RNAforesterOptions::LocalSimilarity) && options.has(RNAforesterOptions::LocalSubopts))
			  cout << "calculate suboptimals within " << suboptPercent << "% of global optimum" << endl << endl;

		// profile search
		if(options.has(RNAforesterOptions::ProfileSearch))
		  {
		    string filename;
		    
		    options.get(RNAforesterOptions::ProfileSearch,filename,string(""));
		    if(filename=="")
		      {
			cerr << "no profile filename" << endl;
			exit(EXIT_FAILURE);
		      }
		    
		    RNAProfileAlignment *rnaProfileAli=new RNAProfileAlignment(filename);
		    alignList.push_back(rnaProfileAli);
		  }

		if(options.has(RNAforesterOptions::ReadFromFile))
		{
			string filename;
			options.get(RNAforesterOptions::ReadFromFile,filename,string(""));
			inputFile=new ifstream(filename.c_str());
			if(inputFile->fail())
			{
				cerr << "cannot open file: \"" << filename << "\"" << endl;
				exit(EXIT_FAILURE);
			}
			inputStream=inputFile;
		}
		else
		{
			inputStream=&cin;
		}


		if(showScale)
		  {
		    if(!options.has(RNAforesterOptions::NoScale) && !options.has(RNAforesterOptions::ReadFromFile))
		      cout << endl << PROMPT << SCALE;
		    showScale=false;
		  }

		for(;;)
		{
			getline(*inputStream,buffer);

			if(inputStream->eof())
			  {			    
			    if(options.has(RNAforesterOptions::Multiple) && !options.has(RNAforesterOptions::ProfileSearch))
				buffer="&";
			    else
			      exit(EXIT_SUCCESS);
			  }

			if(buffer.empty())
				continue;

			// quit if character is @
			if(buffer[0]=='@')
				break;

			// delete '\r' at line end from non unix files
			if(buffer[buffer.size()-1]=='\r')
				buffer.erase(buffer.size()-1);

			// check for name of structure
			if(buffer[0]=='>')
			{
				nameStr=&buffer[1];
				continue;
			}

			// cut after blank
			cutAfterChar(buffer,' ');

			// check for aligning multiple structures
			// if input is read from file the eof has the same meaning as &
			if( buffer[0]=='&')
				multipleAlign=true;
			else
			{
				// check for base string
			  
			  if(RNAFuncs::isRNAString(buffer))
				{
				  baseStr=buffer;
				  // convert to small letters
				  transform(baseStr.begin(),baseStr.end(),baseStr.begin(),ToLower());
				  // t -> u
				  replace(baseStr.begin(),baseStr.end(),'t','u');
				  // delete '.'  and '-' from alignment files
				  remove(baseStr.begin(),baseStr.end(),'.');
				  remove(baseStr.begin(),baseStr.end(),'-');

#ifdef HAVE_LIBRNA  // This features require the ViennaRNA library				  
				  if(options.has(RNAforesterOptions::PredictProfile))
				    {
				      ostringstream ss;		 
				      string constraint;

				      // if there is no name given (> ...) use a counter
				      if(nameStr=="")
					ss << "> " << structure_count;
				      else
					ss << nameStr;	
				      
				      cout << "Predicting structure profile for sequence: " << ss.str() << endl;
				      RNAProfileAlignment *rnaProfileAli=new RNAProfileAlignment(baseStr,ss.str(),constraint,minPairProb);
				      alignList.push_back(rnaProfileAli);
				      makeDotFileInp(*rnaProfileAli,options,structure_count);
				      structure_count++;
				    }
#endif				  

				  //				  continue;
				}
			  else
			    {

				// check for vienna string	
			  if(RNAFuncs::isViennaString(buffer,basePairCount,maxDepth))
			    {

#ifdef HAVE_LIBRNA  // This features require the ViennaRNA library			      
			      // skip structure lines if structures are predicted
			      if(options.has(RNAforesterOptions::PredictProfile))
				{
				  cout << "ignoring structure: " << buffer << endl;
				  continue;
				}
#endif
			      
			     viennaStr=buffer;	       
			    }
			  else   
				{ 
					cerr << "The input sequence is neither an RNA/DNA string nor in vienna format." << endl;
					cerr << "line: " << buffer << endl;
					showScale=true;
					exit(EXIT_FAILURE);
				}

				// add structure to input list
				if(options.has(RNAforesterOptions::Multiple))
				{
					ostringstream ss;		 

					// if there is no name given (> ...) use a counter
					if(nameStr=="")
					  ss << "> " << structure_count;
					else
					  ss << nameStr;	


					RNAProfileAlignment *rnaProfileAli=new RNAProfileAlignment(baseStr,viennaStr,ss.str());
					makeDotFileInp(*rnaProfileAli,options,structure_count);					
					alignList.push_back(rnaProfileAli); 	
				}
				else
				{
				  if(options.has(RNAforesterOptions::TreeEdit))
				    {
				      RNAForestSZ *rnaForestSZ=new RNAForestSZ(baseStr,viennaStr,nameStr);
				      inputListSZ.push_back(rnaForestSZ);
				    }
				  else
				    {
				      RNAForest *rnaForest=new RNAForest(baseStr,viennaStr,nameStr);			
				      nameStr="";
				      makeDotFileInp(*rnaForest,options,structure_count);
				      inputListPW.push_back(rnaForest);
				    }
				}

				structure_count++;      
				showScale=true;
			}
			}

			// ***** multiple alignment
			if((options.has(RNAforesterOptions::Multiple) && multipleAlign) || (options.has(RNAforesterOptions::ProfileSearch) && alignList.size()==2))
			{
				alignMultiple(alignList,score,options);
				multipleAlign=false;
				structure_count=1;
				
				if(options.has(RNAforesterOptions::ProfileSearch))
				  {
				    string filename;
		    
				    options.get(RNAforesterOptions::ProfileSearch,filename,string(""));		    		    
				    RNAProfileAlignment *rnaProfileAli=new RNAProfileAlignment(filename);
				    alignList.push_back(rnaProfileAli);
				  }
				else
				  break;

				//				break;
			}

			// ***** pairwise alignment
			if(inputListPW.size()==2)
			{
			  if(options.has(RNAforesterOptions::GlobalAlignment))
			    alignPairwiseSimple(inputListPW,score,options);
			  else
			    alignPairwise(inputListPW,score,options);
			  break;
			}
			
			if(inputListSZ.size()==2)
			{
			  editPairwise(inputListSZ,score,options);
			  break;
			}

		}

		// free dynamic allocated memory
		deque<RNAForest*>::const_iterator it;
		for(it = inputListPW.begin(); it!=inputListPW.end(); it++)
			delete *it;

		DELETE(inputFile);

		//	  getchar();	// only for testing
		return (0);
	}
	catch(RNAforesterOptions::IncompatibleException e)
	{
		e.showError();
		return(EXIT_FAILURE);
	} 
	catch(RNAforesterOptions::RequiresException e)
	{
		e.showError();
		return(EXIT_FAILURE);
	} 
}



void cutAfterChar(string &s,char c){
	string::size_type pos=s.find(c);
	if(pos!=string::npos)
		s.erase(pos);
}

void alignMultiple(deque<RNAProfileAlignment*> &alignList, Score &score,const RNAforesterOptions &options)
{
	DoubleScoreProfileAlgebraType *alg;
	deque<pair<double,RNAProfileAlignment*> > resultList;
//	double optScore;
	Uint clusterNr=1;
	double minPairProb;

	options.get(RNAforesterOptions::ConsensusMinPairProb,minPairProb,0.5);
	
	// distance or similarity
	if(options.has(RNAforesterOptions::CalculateDistance))
		alg=new DoubleDistProfileAlgebra(score);
	else
		alg=new DoubleSimiProfileAlgebra(score);

	cout << endl;	

	progressiveAlign(alignList,resultList,alg,options);

	cout << endl;
	cout << "*** Results ***" << endl << endl; 
	cout << "Minimum basepair probability for consensus structure (-cmin): " << minPairProb << endl << endl;

	deque<pair<double,RNAProfileAlignment*> >::const_iterator it;
	for(it=resultList.begin();it!=resultList.end();it++)
	  {
	    cout << "RNA Structure Cluster Nr: " << clusterNr << endl;
	    cout << "Score: " << it->first << endl;
	    cout << "Members: " << it->second->getNumStructures() << endl << endl;

	    if(options.has(RNAforesterOptions::FastaOutput))