rnafuncs.cpp

ViennaRNA-1.6.1

      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);

#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);
#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);
  g2_set_line_width(id,0.2);

  // mark 5' end
  g2_string(id,X[0]-20,Y[0],"5'");

  // define colors
  if(options.greyColors)
    {
      ps_color_black=g2_ink(id_PS,0,0,0);
      ps_color_red=g2_ink(id_PS,0,0,0);
      ps_color_blue=g2_ink(id_PS,0.7,0.7,0.7);
      ps_color_green=g2_ink(id_PS,0.4,0.4,0.4);

#ifdef HAVE_LIBGD
      if(options.generatePNG)
	{
	  png_color_black=g2_ink(id_PNG,0,0,0);
	  png_color_red=g2_ink(id_PNG,0,0,0);
	  png_color_blue=g2_ink(id_PNG,0.7,0.7,0.7);
	  png_color_green=g2_ink(id_PNG,0.4,0.4,0.4);
	}

       if(options.generateJPG)
	{
	  jpg_color_black=g2_ink(id_JPG,0,0,0);
	  jpg_color_red=g2_ink(id_JPG,0,0,0);
	  jpg_color_blue=g2_ink(id_JPG,0.7,0.7,0.7);
	  jpg_color_green=g2_ink(id_JPG,0.4,0.4,0.4);  
	}     
#endif
    }
  else
    {
      ps_color_black=g2_ink(id_PS,0,0,0);
      ps_color_red=g2_ink(id_PS,1,0,0);
      ps_color_blue=g2_ink(id_PS,0,0,0.5);
      ps_color_green=g2_ink(id_PS,0,0.5,0);

#ifdef HAVE_LIBGD
      if(options.generatePNG)
	{
	  png_color_black=g2_ink(id_PNG,0,0,0);
	  png_color_red=g2_ink(id_PNG,1,0,0);
	  png_color_blue=g2_ink(id_PNG,0,0,0.5);
	  png_color_green=g2_ink(id_PNG,0,0.5,0);
	}

      if(options.generateJPG)
	{
	  jpg_color_black=g2_ink(id_JPG,0,0,0);
	  jpg_color_red=g2_ink(id_JPG,1,0,0);
	  jpg_color_blue=g2_ink(id_JPG,0,0,0.5);
	  jpg_color_green=g2_ink(id_JPG,0,0.5,0);
	}
#endif
    }  

  // draw sequence
  g2_set_font_size(id,base_fontsize*options.scale);
    for(i=0;i<structure.size();i++)
   {
     isDel = false;
     isIns = false;

     //base
     g2_pen(id,ps_color_black); // match 
     base = "";
     if(seq1[i]!='-')
       {
	 base += seq1[i];
	 basenr_x++;
       }
     else
       {
	 g2_pen(id,ps_color_green); // insertion
	 isIns=true;
       }
     
     if(seq2[i]!='-')
       {
	 base += seq2[i];
	 basenr_y++;
       }
     else
       {
	 g2_pen(id,ps_color_blue); // deletion 
	 isDel=true;
       }
     
     if(base.size()==2 && base[0]!=base[1])
       g2_pen(id,ps_color_red); // mismatch  
     else
       base = base[0];          // show duplicate base only once
     
     xpos=X[i]-base.length()*(base_fontsize*options.scale)/2;
     ypos=Y[i]-4;
     g2_string(id,xpos,ypos,(char*)base.c_str());

     if(!options.hideBaseNumbers)
       {
	 // draw base number
	 if(!isIns && basenr_x % options.baseNumInterval == 0)
	   {	 
	     g2_pen(id,ps_color_blue);
	     sprintf(buf,"%d",basenr_x);
	     g2_string(id,xpos-20,ypos,buf);
	   }
	 if(!isDel && basenr_y % options.baseNumInterval == 0)
	   {	 
	     g2_pen(id,ps_color_green);
	     sprintf(buf,"%d",basenr_y);
	     g2_string(id,xpos+20,ypos,buf);
	   }
       }
     
     // connection to next base
     if(i<structure.size()-1)
       {
	 // if the connected bases are only in one of the structures
	 // use the appropriate color
	 if(seq1[i]=='-' && seq1[i+1]=='-')
	   g2_pen(id,ps_color_green);
	 else
	   if(seq2[i]=='-' && seq2[i+1]=='-')
	     g2_pen(id,ps_color_blue);
	   else	      
	     g2_pen(id,ps_color_black);

	 g2_line(id,X[i],Y[i],X[i+1],Y[i+1]);

	 // draw circles at line endpoints
	 if(seq1[i]=='-')
	   g2_pen(id,ps_color_green);
	 else
	   if(seq2[i]=='-')
	     g2_pen(id,ps_color_blue);
	   else	      
	     g2_pen(id,ps_color_black); 
	 g2_filled_circle(id,X[i],Y[i],0.7*options.scale);       // circles are drawn twice, but thats ok ...

	 if(seq1[i+1]=='-')
	   g2_pen(id,ps_color_green);
	 else
	   if(seq2[i+1]=='-')
	     g2_pen(id,ps_color_blue);
	   else	      
	     g2_pen(id,ps_color_black);
	 g2_filled_circle(id,X[i+1],Y[i+1],0.7*options.scale);
       }
   }
   
    // draw pairings
    // !!! pair_table indexing begins at 1 !!!
    for(i=0;i<structure.size();i++)
      {
	if((unsigned short)pair_table[i+1]>i+1 && (unsigned short)alt_pair_table[i+1]>i+1)
	  {	    	    
	    // pairs in both structures
	    g2_pen(id,ps_color_black);	   	    
	    g2_line(id,X[i],Y[i],X[pair_table[i+1]-1],Y[pair_table[i+1]-1]);
	  }
	else
	  {
	    if((unsigned short)pair_table[i+1]>i+1 && (unsigned short)alt_pair_table[i+1]<=i+1)
	      {	    	    
		// pairs only in first structure
		if(atX)
		  g2_pen(id,ps_color_blue);
		else
		  g2_pen(id,ps_color_green);

		g2_line(id,X[i],Y[i],X[pair_table[i+1]-1],Y[pair_table[i+1]-1]);
	      }
	    else
	      {
		if((unsigned short)pair_table[i+1]<=i+1 && (unsigned short)alt_pair_table[i+1]>i+1)
		  {	    	    
		    // pairs only in second structure
		    if(atX)
		      g2_pen(id,ps_color_green);	   	    
		    else
		      g2_pen(id,ps_color_blue);

		    double dashes=2.0;
		    g2_set_dash(id,1,&dashes);
		    g2_line(id,X[i],Y[i],X[alt_pair_table[i+1]-1],Y[alt_pair_table[i+1]-1]);
		    dashes=0;
		    g2_set_dash(id,1,&dashes);
		  }
	      }
	  }
      }


    // draw structure names
    // x-at-y or y-at-x
    if(atX)
      {
	g2_pen(id,ps_color_green);
	g2_string(id,min_X,max_Y-10,(char*)strname2.c_str());
	g2_pen(id,ps_color_black);
	g2_string(id,min_X,max_Y-20,"at");
	g2_pen(id,ps_color_blue);
	g2_string(id,min_X,max_Y-30,(char*)strname1.c_str());
      }
    else
      {
	g2_pen(id,ps_color_blue);
	g2_string(id,min_X,max_Y-10,(char*)strname1.c_str());
	g2_pen(id,ps_color_black);
	g2_string(id,min_X,max_Y-20,"at");
	g2_pen(id,ps_color_green);
	g2_string(id,min_X,max_Y-30,(char*)strname2.c_str());
      }	   

    g2_flush(id);
    g2_close(id);

    free(pair_table);
    free(alt_pair_table);
    DELETE(X);
    DELETE(Y);
}
#endif

#ifdef HAVE_LIBRNA
void RNAFuncs::generateRNAAlignmentXML(const string &structure, const string &altStructure, const string &seq1, const string &seq2, const string &strname1, const string &strname2, ostream &s)
{
  string base;
  Uint i;
  float *X, *Y;
  short *pair_table,*alt_pair_table;
  Uint basenr_x=1, basenr_y=1;
  string filename;

  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());

  // calculate coordinates
  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;
    
  // generate XML
  s << "<alignment xname=\"" << strname1 << "\" yname=\"" << strname2 << "\">" << endl;

  // seqalignment
  s << "  <seqalignment>" << endl;
  for(i=0;i<seq1.length();i++)
    {
      s << "    <base id=\"" << i+1 << "\" xbasenr=\"" << basenr_x << "\" ybasenr=\"" << basenr_y << "\" xbase=\"" << seq1[i] <<"\" ybase=\"" << seq2[i] << "\" />" << endl;
      if(seq1[i]!='-')
	basenr_x++;

      if(seq2[i]!='-')
	basenr_y++;      
    }
  s << "  </seqalignment>" << endl;

  // pairs
  s << "  <pairs>" << endl;
  for(i=0;i<structure.size();i++)
    {
      // both
      if((unsigned short)pair_table[i+1]>i+1 && (unsigned short)alt_pair_table[i+1]>i+1) 
	{
	  s << "    <pair drawbaseid1=\"" << i+1 << "\" drawbaseid2=\"" << pair_table[i+1] << "\"/>" << endl;
	}
      else
	{
	  if((unsigned short)pair_table[i+1]>i+1)
	    {
	      s << "    <pair drawbaseid1=\"" << i+1 << "\" drawbaseid2=\"" << pair_table[i+1] << "\" structid=\"1\"/>" << endl;
	    }
	  else
	    {
	    if((unsigned short)alt_pair_table[i+1]>i+1)
	      s << "    <pair drawbaseid1=\"" << i+1 << "\" drawbaseid2=\"" << alt_pair_table[i+1] << "\" structid=\"2\"/>" << endl;
	    }
	}
    }
  s << "  </pairs>" << endl;

  // coordinates
  // x structure
  s << "  <structure id=\"1\">" << endl;
  s << "    <drawbases>" << endl;
  for(i=0;i<seq1.length();i++)
    {
      s << "      <drawbase id=\"" << i+1 << "\" xcoor=\"" << X[i] << "\" ycoor=\"" << -Y[i] << "\"/>" << endl;
    }
  s << "    </drawbases>" << endl;
  s << "  </structure>" << endl;
  
  DELETE(X);
  DELETE(Y);
  X = new float[altStructure.size()];
  Y = new float[altStructure.size()];

  // y structure
  i = naview_xy_coordinates(alt_pair_table, X, Y);
  if(i!=structure.size())
    cerr << "strange things happening in squigglePlot ..." << endl;

  s << "  <structure id=\"2\">" << endl;
  s << "    <drawbases>" << endl;
  for(i=0;i<seq1.length();i++)
    {
      s << "      <drawbase id=\"" << i+1 << "\" xcoor=\"" << X[i] << "\" ycoor=\"" << -Y[i] << "\"/>" << endl;
    }
  s << "    </drawbases>" << endl;
  s << "  </structure>" << endl;
  s << "</alignment>" << endl;

  free(pair_table);
  free(alt_pair_table);
  DELETE(X);
  DELETE(Y);
}
#endif

void RNAFuncs::printAli(const string &name1, const string &name2, const string &seq1, const string &seq2, const string &str1, const string &str2)
{
  Uint i;
  string info;

  for(i=0;i<seq1.length();i++)		 
    info += seq1[i]==seq2[i] ? "*" : " ";		  		

  for(i=0;i<seq1.length();i+=55)
    {
      cout << setw(20) << setfill(' ') << left << name1.substr(0,20) << setw(5) << " " << seq1.substr(i,55) << endl;
      cout << setw(20) << setfill(' ') << left << name2.substr(0,20) << setw(5) << " " << seq2.substr(i,55) << endl;
      cout << setw(25) << " " << info.substr(i,55) << endl;
    }

  cout << endl;

  info="";
  for(i=0;i<str1.length();i++)		 
    info += str1[i]==str2[i] ? "*" : " ";

  for(i=0;i<str1.length();i+=55)
    {
      cout << setw(20) << setfill(' ') << left << name1.substr(0,20) << setw(5) << " " << str1.substr(i,55) << endl;
      cout << setw(20) << setfill(' ') << left << name2.substr(0,20) << setw(5) << " " << str2.substr(i,55) << endl;
      cout << setw(25) << " " << info.substr(i,55) << endl;
    }

  cout << endl;

}

Uint RNAFuncs::treeSize(const string &viennaStr)
{
   Ulong basePairCount=0,maxDepth=0;

   RNAFuncs::isViennaString(viennaStr,basePairCount,maxDepth);
   // the size of the forests is determined by the number of bases and basepairs	
   return viennaStr.length() + basePairCount;  
};