rnafuncs.cpp

ViennaRNA-1.6.1

#include <algorithm>
#include <cctype>
#include <iomanip>
#include <iostream>
#include <string>

#include "debug.h"
#include "misc.h"
#include "rnafuncs.h"
#include "utils.h"

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

#ifndef HAVE_LIBRNA
#undef HAVE_LIBG2
#endif

#include <assert.h>

#ifdef HAVE_LIBG2
#include <g2.h>
#include <g2_PS.h>
#include <g2_FIG.h>

#ifdef HAVE_LIBGD
#include <g2_gd.h>
#endif
#endif

#define COLOR_RED       0
#define COLOR_GREEN     1
#define COLOR_BLUE      2
#define COLOR_YELLOW    3
#define COLOR_MAGENTA   4
#define COLOR_TURKIS    5
#define NUM_COLORS      6

#define COLOR_DEF_RED       1,0,0
#define COLOR_DEF_GREEN     0,1,0
#define COLOR_DEF_BLUE      0,0,1
#define COLOR_DEF_YELLOW    1,1,0
#define COLOR_DEF_MAGENTA   0,1,1
#define COLOR_DEF_TURKIS    1,0,1

bool RNAFuncs::isRNAString(const string &str)
{	
  string::const_iterator it;
  for(it=str.begin(); it!=str.end(); it++)
    {
      switch(tolower(*it))
	{
	case 'a':
	case 'c':
	case 'g':
	case 'u':
	case 't':
	case '.':
	case '-':
	  break;
	default:
	  return false;
	}
    }

  return true;
}


bool RNAFuncs::isViennaString(const string &str, Ulong &basePairCount, Ulong &maxDepth)
{
  long brackets=0;
  maxDepth=0;
  basePairCount=0;       

  string::const_iterator it;
  for(it=str.begin(); it!=str.end(); it++)
    {
      switch(*it)
	{
	case '.':
	  break;
	case '(':
	  brackets++;
	  basePairCount++;
	  if(brackets>0)	
	    maxDepth=max(static_cast<Ulong>(maxDepth),static_cast<Ulong>(brackets));
	  break;
	case ')':
	  if(brackets)
	    {
	      brackets--;
	      break;
	    }
	default:
	  return false;
	}
    }

  // equal nr of opening and closing brackets ?
  if(brackets)
    return false;

  return true;
}

#ifdef HAVE_LIBG2  // This features require the g2 library
void RNAFuncs::drawRNAStructure(const string &seq, const string &structure, const string &filename_prefix, const string &structname, const list<pair<Uint,Uint> > &regions, const SquigglePlotOptions &options)
{
  const double base_fontsize=12;

  float *X, *Y,min_X=0,max_X=0,min_Y=0,max_Y=0;
  Uint i;
  short *pair_table;
  int id_PS,id_FIG=0,id;
  int color_black, *colors;
  double xpos,ypos;
  char buf[10];
  string filename;
  double *points;
  int numPoints;

#ifdef HAVE_LIBGD
  int id_PNG=0,id_JPG=0;
#endif

  X = new float[structure.size()];
  Y = new float[structure.size()];
  points=new double[2*structure.size()];
  colors=new int[NUM_COLORS];

  assert(seq.size() == structure.size());

  pair_table = make_pair_table(structure.c_str());
  i = naview_xy_coordinates(pair_table, X, Y);
  if(i!=structure.size())
    cerr << "strange things happening in squigglePlot ..." << endl;
    
  // scale image
  for(i=0;i<structure.size();i++)
    {
      X[i]*=static_cast<float>(options.scale);
      Y[i]*=static_cast<float>(options.scale);
    }  

  // calculate image dimensions
  for(i=0;i<structure.size();i++)
    {
      min_X=min(min_X,X[i]);
      max_X=max(max_X,X[i]);
      min_Y=min(min_Y,Y[i]);
      max_Y=max(max_Y,Y[i]);
    }

  // add a border to image size
  min_X-=10;
  max_X+=10;
  min_Y-=10;
  max_Y+=10;

  //id_PS  = g2_open_PS("ali.ps", g2_A4, g2_PS_port);
  filename=filename_prefix + ".ps";
  id_PS  = g2_open_EPSF((char*)filename.c_str());
  g2_set_coordinate_system(id_PS,-min_X,-min_Y,1,1);

  if(options.generateFIG)
    {
      filename=filename_prefix + ".fig";
      id_FIG=g2_open_FIG((char*)filename.c_str());
      g2_set_coordinate_system(id_FIG,-min_X,-min_Y,1,1);
    }
#ifdef HAVE_LIBGD
  if(options.generatePNG)
    {
      filename=filename_prefix + ".png";
      id_PNG=g2_open_gd((char*)filename.c_str(),(int)(max_X-min_X),(int)(max_Y-min_Y),g2_gd_png);
      g2_set_coordinate_system(id_PNG,-min_X,-min_Y,1,1);
    }
   if(options.generateJPG)
     {
       filename=filename_prefix + ".jpg"; 
       id_JPG=g2_open_gd((char*)filename.c_str(),(int)(max_X-min_X),(int)(max_Y-min_Y),g2_gd_jpeg);
       g2_set_coordinate_system(id_PS,-min_X,-min_Y,1,1);
     }
#endif

  id     = g2_open_vd();
  g2_attach(id,id_PS);

  if(options.generateFIG)
    {
      g2_attach(id,id_FIG);
    }

#ifdef HAVE_LIBGD
  if(options.generatePNG)
    g2_attach(id,id_PNG);
  if(options.generateJPG)
    g2_attach(id,id_JPG);
#endif

  // cout << "min_X: " << min_X <<",max_X: " << max_X << ",min_Y: " << min_Y << "max_Y: " << max_Y << endl; 
  //  g2_set_coordinate_system(id_PS,595/2.0,842/2.0,0.5,0.5);

  // define colors
  if(options.greyColors)
    {
      color_black=g2_ink(id_PS,0,0,0);
      for(i=0;i<NUM_COLORS;i++)
	colors[i]=g2_ink(id_PS,0,0,0);

      if(options.generateFIG)
	{
	  color_black=g2_ink(id_FIG,0,0,0);
	  for(i=0;i<NUM_COLORS;i++)
	    colors[i]=g2_ink(id_FIG,0,0,0);
	}
      
#ifdef HAVE_LIBGD
      if(options.generatePNG)
	{
	  color_black=g2_ink(id_PNG,0,0,0);
	  for(i=0;i<NUM_COLORS;i++)
	    colors[i]=g2_ink(id_PNG,0,0,0);	  
	}

      if(options.generateJPG)
	{
	  color_black=g2_ink(id_JPG,0,0,0);
	  for(i=0;i<NUM_COLORS;i++)
	    colors[i]=g2_ink(id_JPG,0,0,0);	  
	}     
#endif
    }
  else
    {
      color_black=g2_ink(id,0,0,0);
      colors[COLOR_RED]=g2_ink(id_PS,COLOR_DEF_RED);
      colors[COLOR_GREEN]=g2_ink(id_PS,COLOR_DEF_GREEN);
      colors[COLOR_BLUE]=g2_ink(id_PS,COLOR_DEF_BLUE);
      colors[COLOR_YELLOW]=g2_ink(id_PS,COLOR_DEF_YELLOW);
      colors[COLOR_MAGENTA]=g2_ink(id_PS,COLOR_DEF_MAGENTA);
      colors[COLOR_TURKIS]=g2_ink(id_PS,COLOR_DEF_TURKIS);

      if(options.generateFIG)
	{
	  color_black=g2_ink(id_FIG,0,0,0);
	  colors[COLOR_RED]=g2_ink(id_FIG,COLOR_DEF_RED);
	  colors[COLOR_GREEN]=g2_ink(id_FIG,COLOR_DEF_GREEN);
	  colors[COLOR_BLUE]=g2_ink(id_FIG,COLOR_DEF_BLUE);
	  colors[COLOR_YELLOW]=g2_ink(id_FIG,COLOR_DEF_YELLOW);
	  colors[COLOR_MAGENTA]=g2_ink(id_FIG,COLOR_DEF_MAGENTA);
	  colors[COLOR_TURKIS]=g2_ink(id_FIG,COLOR_DEF_TURKIS);
	}
      
#ifdef HAVE_LIBGD
      if(options.generatePNG)
	{
	  color_black=g2_ink(id_PNG,0,0,0);
	  colors[COLOR_RED]=g2_ink(id_PNG,COLOR_DEF_RED);
	  colors[COLOR_GREEN]=g2_ink(id_PNG,COLOR_DEF_GREEN);
	  colors[COLOR_BLUE]=g2_ink(id_PNG,COLOR_DEF_BLUE);
	  colors[COLOR_YELLOW]=g2_ink(id_PNG,COLOR_DEF_YELLOW);
	  colors[COLOR_MAGENTA]=g2_ink(id_PNG,COLOR_DEF_MAGENTA);
	  colors[COLOR_TURKIS]=g2_ink(id_PNG,COLOR_DEF_TURKIS);
	}

      if(options.generateJPG)
	{
	  color_black=g2_ink(id_JPG,0,0,0);
	  colors[COLOR_RED]=g2_ink(id_JPG,COLOR_DEF_RED);
	  colors[COLOR_GREEN]=g2_ink(id_JPG,COLOR_DEF_GREEN);
	  colors[COLOR_BLUE]=g2_ink(id_JPG,COLOR_DEF_BLUE);
	  colors[COLOR_YELLOW]=g2_ink(id_JPG,COLOR_DEF_YELLOW);
	  colors[COLOR_MAGENTA]=g2_ink(id_JPG,COLOR_DEF_MAGENTA);
	  colors[COLOR_TURKIS]=g2_ink(id_JPG,COLOR_DEF_TURKIS);
	}
#endif
    }
  
  for(i=0;i<structure.size();i++)
    {
      // connection to next base
      if(i<structure.size()-1)
	{
	 // if the connected bases are only in one of the structures
	  // use the appropriate color
	  g2_line(id,X[i],Y[i],X[i+1],Y[i+1]);
	  
	  // draw circles at line endpoints
	  g2_filled_circle(id,X[i],Y[i],0.7*options.scale);       // circles are drawn twice, but thats ok ...
	}
    }

  // draw pairings
  // !!! pair_table indexing begins at 1 !!!
  for(i=0;i<structure.size();i++)
    {
	if((unsigned short)pair_table[i+1]>i+1)
	  {	    	    
	    // pairs in both structures
	    g2_line(id,X[i],Y[i],X[pair_table[i+1]-1],Y[pair_table[i+1]-1]);
	  }
    }

  // draw regions
  g2_set_line_width(id,0.4);
  list<pair<Uint,Uint> >::const_iterator it;  
  Uint regionNr=0;

  for(it=regions.begin();it!=regions.end();it++)
    {
      double center_x=0,center_y=0;

      // fill coordinate list
	for(i=0;i < it->second;i++)
	  {
	    points[2*i]=X[it->first+i];
	    points[2*i+1]=Y[it->first+i];

	    center_x+=X[it->first+i];
	    center_y+=Y[it->first+i];
	  }
	numPoints=it->second-1;
	center_x/=numPoints;   // center of gravity
	center_y/=numPoints;
	
	g2_pen(id,colors[regionNr % NUM_COLORS]);
	g2_poly_line(id,numPoints,points);
	sprintf(buf,"%d",regionNr+1);
	g2_string(id,center_x,center_y,buf);   // draw region number

	regionNr++;
    }
  
  // mark 5' end

  g2_pen(id,color_black);
  g2_string(id,X[0]-20,Y[0],"5'");

  g2_set_font_size(id,base_fontsize*options.scale);
  g2_set_line_width(id,0.2);
  
  // draw sequence
  for(i=0;i<structure.size();i++)
    {
     g2_pen(id,color_black); // match 
     xpos=X[i]-(base_fontsize*options.scale)/2;
     ypos=Y[i]-4;
     sprintf(buf,"%c",seq[i]);
     g2_string(id,xpos,ypos,buf);

     /*
     if(!options.hideBaseNumbers)
       {
	 // draw base number
	 if(basenr % options.baseNumInterval == 0)
	   {	 
	     sprintf(buf,"%d",basenr);
	     g2_string(id,xpos-20,ypos,buf);
	   }
	   }*/     
    }

    // draw structure name
    g2_string(id,min_X,max_Y-10,(char*)structname.c_str());
        
    g2_flush(id);
    g2_close(id);

    free(pair_table);
    DELETE(X);
    DELETE(Y);
    DELETE(points);
    DELETE(colors);
}
#endif

#ifdef HAVE_LIBG2  // This features require the g2 library
void RNAFuncs::drawRNAAlignment(const string &structure, const string &altStructure, const string &seq1, const string &seq2, const string &strname1, const string &strname2, const string &filename_prefix, const bool atX, const SquigglePlotOptions &options)
{
  const double base_fontsize=12;

  string base,structname;
  float *X, *Y,min_X=0,max_X=0,min_Y=0,max_Y=0;
  Uint i;
  short *pair_table,*alt_pair_table;
  int id_PS,id;
  int ps_color_black,ps_color_red,ps_color_blue,ps_color_green;
  Uint basenr_x=0, basenr_y=0;
  double xpos,ypos;
  char buf[10];
  bool isDel,isIns;
  string filename;
#ifdef HAVE_LIBGD
  int id_PNG=0,id_JPG=0;
  int png_color_black,png_color_red,png_color_blue,png_color_green;
  int jpg_color_black,jpg_color_red,jpg_color_blue,jpg_color_green;
#endif

  X = new float[structure.size()];
  Y = new float[structure.size()];

  assert(seq1.size() == seq2.size());
  assert(seq1.size() == structure.size());
  assert(seq1.size() == altStructure.size());

  pair_table = make_pair_table(structure.c_str());
  alt_pair_table = make_pair_table(altStructure.c_str());
  i = naview_xy_coordinates(pair_table, X, Y);
  if(i!=structure.size())
    cerr << "strange things happening in squigglePlot ..." << endl;
    
  // scale image
  for(i=0;i<structure.size();i++)
    {
      X[i]*=static_cast<float>(options.scale);