dna2.cc

unix或linux下的DNA分析软件源码 其功能如下 1. Edit up to 256 peptide or DNA sequences simultaneously. 2. Transla

//--------------------------------------------------------------------------//
// dna2.cc                                                                  //
// Latest revision: 05-06-2000                                              //
// Copyright (C) 2000 by Thomas J. Nelson                                   //
// All rights reserved.                                                     //
//--------------------------------------------------------------------------//

#include "dna.h"

extern Globals g;
//----Module-level globals ----------------------------------//  
int ylabel[30] = { 220, 235, 250, 265, 280, 295, 310, 325, 340, 355, 370,
                   385, 400, 415, 430, 445, 460, 475, 490 };    

/* IMPORTANT */
/* NOTE: If 'text' is a DNA sequence, any changes in 'text' must be 
   copied to 'dtext'. dtext must be kept current in case user wants
   to translate to a different reading frame. */

//--------------------------------------------------------------------------//
// analyze - compare two sequences                                          //
//--------------------------------------------------------------------------//
void analyze(void)
{
  int n1, n2;
  g.editing = 0;
  n1 = g.doc[g.matrix_set1].n;
  n2 = g.doc[g.matrix_set2].n;
  int b = (7+g.bitsperpixel)/8;
  if(g.image_1d != NULL) free(g.image_1d);
  g.image_1d = (uchar*)malloc(n1*n2*sizeof(uchar)*b);
  if(g.image_1d==0) { fprintf(stderr,"%s for image_1d\n",g.nomemory); return; }
  g.image = make_2d_alias(g.image_1d, n1*b, n2);
  g.image_ximage = createximage(n1, n2, g.bitsperpixel, g.image_1d);
  compare(g.doc[g.matrix_set1].text, g.doc[g.matrix_set2].text, n1, n2, g.image);
  redrawscreen(0, 0, n1, n2);
}



//--------------------------------------------------------------------------//
// align                                                                    //
//--------------------------------------------------------------------------//
double align(void)
{
  int bestoffset=0, n1, n2, offset;
  n1 = g.doc[0].n;
  n2 = g.doc[1].n;
  double score[2*n1], bestscore=0;
  for(offset=-n1; offset<n1; offset++)
  {   score[offset+n1] = alignment(g.doc[0].text, g.doc[1].text, offset);
      if(score[offset+n1] > bestscore) 
      {  bestscore = score[offset+n1];
         bestoffset = offset;
      }
  } 

  // Make sure offset is positive

  if(bestoffset<0){ g.doc[0].offset = 0; g.doc[1].offset = -bestoffset; }
  else            { g.doc[0].offset = bestoffset; g.doc[1].offset = 0; }
  g.doc[0].pos = 0;
  g.doc[1].pos = 0;
  g.doc[0].row = 0;
  g.doc[1].row = 0;
  XClearWindow(g.display, g.draw_window);
  redrawscreen(0, 0, g.xsize, g.ysize);
  return bestscore;
}


//--------------------------------------------------------------------------//
// alignment  - calculate % identity of 2 sequences                         //
// Sequences must contain only valid bases.                                 //
//--------------------------------------------------------------------------//
double alignment(char *s1, char *s2, int offset)
{     
   if(s1==NULL || s2==NULL) return 0;
   double f = 0.0;
   int k, matches=0, n1, n2, total=0;
   int start, end;

   n1 = (int)strlen(s1) + offset;
   n2 = (int)strlen(s2);
   start = max(offset,0);
   if(g.alignsize==0) end = min(n1,n2);
   else               end = min(start+g.alignsize, min(n1,n2));

   for(k=start;k<end;k++)
       if(k+offset>=0 && k>=0) if(s1[k+offset]==s2[k]) matches++;
   total = min((int)strlen(s1), (int)strlen(s2));
   if(total) f = ((double)matches) / ((double)(total));
   return f;
}



//--------------------------------------------------------------------------//
// compare - compare 2 DNA sequences                                        //
// The strings must be uppercase and contain only valid codes.              //
// Puts the result in a 2D array 'image'.                                   //
//--------------------------------------------------------------------------//
void compare(char *s1, char *s2, int n1, int n2, uchar **image)
{
  int bpp, color, b;
  int i,j;
  bpp = g.bitsperpixel;
  b = (bpp+7)/8;

  int k=0;
  for(j=k+0;j<k+n2;j++)
  for(i=0;i<n1;i++)
  {  if(s1[i]==s2[j] && s1[i]!='N' && s1[i]!='X') color=g.fcolor; 
           else color=g.bcolor;
     putpixelbytes(image[j-k]+i*b,color,bpp);
  }
}



//--------------------------------------------------------------------------//
// printmatches - print no. of matches at left                              //
//--------------------------------------------------------------------------//
void printmatches(int matches, int mismatches)
{  
  int k;
  double ratio = 1;
  char string[256];
  int total = matches + mismatches;
  int y = 12, len;
  Window win = XtWindow(g.match_area);
  XSetForeground(g.display, g.gc, g.fcolor);
  XSetBackground(g.display, g.gc, g.bcolor);
  if(total>0) ratio = (double)matches/(double)total;  
  XDrawImageString(g.display, win, g.gc, 10, y, "Matches:",8);
  sprintf(string,"%d / %d ",matches,total);
  len = (int)strlen(string);
  for(k=len; k<19; k++) string[k]=' ';
  string[19] = 0;
  XDrawImageString(g.display, win, g.gc, 10, y+12, string, 19);
  gcvt(ratio, g.signif, string);   
  len = (int)strlen(string);
  for(k=len; k<19; k++) string[k]=' ';
  string[19] = 0;
  XDrawImageString(g.display, win, g.gc, 10, y+24, string, 19);
}


//--------------------------------------------------------------------------//
// readsequence                                                             //
//--------------------------------------------------------------------------//
void readsequence(void)
{
  int c, n, k, cno=0;
  FILE *fp;
  int fsize;
  char tempstring[1024];
  char junk[1024];
  static char filename[FILENAMELENGTH] = "none";

  //// Set to next sequence no.
  for(k=0; k<MAXSEQUENCES; k++) if(g.doc[k].n==0){ g.cs = k; break; }

  //// Sequence no. is not known until after getfilename()
  strcpy(filename, getfilename(filename));
  if((fp=fopen(filename,"rb"))==NULL) 
  {   message("Can't open file");
     return;
  }
  g.doc[g.cs].ncomments = 0; 

  //// Erase any previous sequence & comments
  init_editor(g.cs);

  // Allocate extra space in case it gets edited
  fsize = filesize(filename);
  fsize = max(65536, 2 * fsize + 10);

  g.doc[g.cs].text  = new char[fsize];  
  g.doc[g.cs].otext = new char[fsize];  
  g.doc[g.cs].dtext = new char[fsize];  
  memset(g.doc[g.cs].text,  0, fsize);
  memset(g.doc[g.cs].otext, 0, fsize);
  memset(g.doc[g.cs].dtext, 0, fsize);
  for(k=0;k<MAXCOMMENTS;k++) g.doc[g.cs].comment[k] = new char[1024];

  if(g.doc[g.cs].text==NULL)
  {  message("Sequence file is too big", ERROR);
     return;
  } 

  n = 0;
  while(!feof(fp)) 
  {  c = toupper(fgetc(fp));
     if(c=='#')
     {   if(cno<MAXCOMMENTS) 
         { 
            fgets(g.doc[g.cs].comment[cno], 1023, fp);  // comment
            cno++;
            g.doc[g.cs].ncomments = cno; 
         }else  
            fgets(junk, 1023, fp);  // comment
     }
     if(valid_base(c, g.sequence_mode)) g.doc[g.cs].text[n++] = c;
  }
  fclose(fp);
  if(g.sequence_mode == PROTEIN)
      sprintf(tempstring,"Read %d amino acids\nfrom file %s\nin sequence no.%d",n,filename,g.cs+1);
  else
      sprintf(tempstring,"Read %d bases\nfrom file %s\nin sequence no.%d",n,filename,g.cs+1);

  message(tempstring);
  strcpy(g.doc[g.cs].filename, filename);
  g.doc[g.cs].type = g.sequence_mode;
  g.doc[g.cs].otype = g.sequence_mode;
  g.doc[g.cs].n = n;
  g.doc[g.cs].on = n;
  g.doc[g.cs].frame = 0;
  memcpy(g.doc[g.cs].otext, g.doc[g.cs].text, n);
  memcpy(g.doc[g.cs].dtext, g.doc[g.cs].text, n);

  g.noofsequences = 0;
  for(k=0; k<MAXSEQUENCES; k++) if(g.doc[k].n) g.noofsequences++;
  XClearWindow(g.display, g.draw_window);
  redrawscreen(0, 0, g.xsize, g.ysize);
  g.touched = 1;
}


//--------------------------------------------------------------------------//
//  init_editor  - initialize values in doc                                 //
//  Call after init_xlib. This should be replaced by a constructor.         //
//--------------------------------------------------------------------------//
void init_editor(int k)
{ 
  int j, direction, ascent, descent;
  XCharStruct overall;

  g.noofsequences = min(g.noofsequences, MAXSEQUENCES);
  g.doc[k].row = 0;
  g.doc[k].pos = 0;
  g.doc[k].text_direction = HORIZONTAL;
  g.doc[k].fontstruct = g.fontstruct;
  XTextExtents(g.doc[k].fontstruct, "Ejyk", 4, &direction, &ascent,
         &descent, &overall);
  g.doc[k].charheight = ascent + descent; 
  g.doc[k].spacing = g.doc[k].charheight;     
  if(g.doc[k].spacing==0) g.doc[k].spacing = 16;

  if(g.doc[k].text != NULL) delete[] g.doc[k].text;
  g.doc[k].text = NULL;

  if(g.doc[k].otext != NULL) delete[] g.doc[k].otext;
  g.doc[k].otext = NULL;

  if(g.doc[k].dtext != NULL) delete[] g.doc[k].dtext;
  g.doc[k].dtext = NULL;

  for(j=0; j<MAXCOMMENTS; j++) 
  {  if(g.doc[k].comment[j] != NULL)
     {   delete[] g.doc[k].comment[j]; 
         g.doc[k].comment[j]=NULL; 
     }
  }
  g.doc[k].ncomments = 0; 

  g.doc[k].touched = 0;
  g.doc[k].marked = 0;
  g.doc[k].markstart = 0;
  g.doc[k].markend = 0;
  if(g.doc[k].filename==NULL)
  {  g.doc[k].filename = new char[FILENAMELENGTH];
     g.doc[k].filename[0] = 0;
  }
  g.doc[k].charwidth = XTextWidth(g.doc[k].fontstruct, "m", 1);
  g.doc[k].topline = 0;
  g.doc[k].linestodisplay = g.main_ysize / g.doc[k].spacing;
  g.doc[k].n = 0;
  g.doc[k].offset = 0;
  g.doc[k].frame = 0;
}


//--------------------------------------------------------------------------//
// move_cursor - move cursor of document k                                  //
//--------------------------------------------------------------------------//
void move_cursor(int k, int row, int pos)
{
  const int ytopmargin = 10;
  static int ocursorx[MAXSEQUENCES], ocursory[MAXSEQUENCES];
  int charsperline, cursorx, cursory, ll;
  charsperline = g.xsize / g.doc[k].charwidth; 
  int dy = g.doc[k].charheight + LINESPACING;

  // Fix cursor position to account for offset

  pos = g.doc[k].pos - g.doc[k].offset;
  row = g.doc[k].row;
  ll = row * charsperline;
  if(pos - ll >= charsperline) { row++; pos -= charsperline; }
  if(pos - ll < 0) { row--; pos += charsperline; }

  cursorx = XTextWidth(g.fontstruct, g.doc[k].text, pos%charsperline);
  cursory = ytopmargin + row * totalspacing() + k * dy;

  XSetForeground(g.display, g.gc, g.bcolor);

  // Erase the old cursor
  XDrawLine(g.display, g.draw_window, g.gc, ocursorx[k], ocursory[k], 
  ocursorx[k] + g.doc[k].charwidth, ocursory[k]);
  XSetForeground(g.display, g.gc, g.fcolor);

  // Draw the new cursor
  XDrawLine(g.display, g.draw_window, g.gc, cursorx, cursory, 
            cursorx + g.doc[k].charwidth, cursory);
  ocursorx[k] = cursorx;
  ocursory[k] = cursory;
}


//--------------------------------------------------------------------------//
// totalspacing - calculate total no. pixels between each sequence          //
//--------------------------------------------------------------------------//
int totalspacing(void)
{
   int k, total=0;  
   int groupspacing = g.doc[0].charheight; // pixels between groups of sequences
   for(k=0; k<g.noofsequences;k++) total += g.doc[k].spacing;
   total += groupspacing; 
   return total;
}




//--------------------------------------------------------------------------//
// redrawscreen                                                             //
// redraw rectangular region of screen                                      //
//--------------------------------------------------------------------------//
int redrawscreen(int xx1, int yy1, int xx2, int yy2)
{
   static int oediting=-1;
   int xx[MAXSEQUENCES], yy[MAXSEQUENCES];
   int charsperline=40, color, charcount=0, dy, dytot, j, k,
      linecount, len, matches=0, mismatches=0, pos, pos2, ascent, descent,
      start, width=0, x, x2, ytopmargin, y2, xshift, xlabelsize, ylabelsize;
   char *s1=NULL, *s2=NULL;
   Window win = XtWindow(g.draw_widget);
   char tempstring[FILENAMELENGTH];

   XSetForeground(g.display, g.gc, g.fcolor);