jpayne@69: //------------------------------------------------------------------------------
jpayne@69: //   Programmer: Adam M Phillippy, The Institute for Genomic Research
jpayne@69: //         File: show-aligns.cc
jpayne@69: //         Date: 10 / 18 / 2002
jpayne@69: //
jpayne@69: //        Usage: show-aligns [options] <deltafile>
jpayne@69: //               Try 'show-aligns -h' for more information
jpayne@69: //
jpayne@69: //  Description: For use in conjunction with the MUMmer package.
jpayne@69: //              "show-aligns" displays human readable information from the
jpayne@69: //             .delta output of the "nucmer" and "promer" programs. Outputs
jpayne@69: //            pairwise alignments to stdout. Works for both nucleotide and
jpayne@69: //           amino-acid alignments.
jpayne@69: //
jpayne@69: //------------------------------------------------------------------------------
jpayne@69: 
jpayne@69: #include "delta.hh"
jpayne@69: #include "tigrinc.hh"
jpayne@69: #include "translate.hh"
jpayne@69: #include "sw_alignscore.hh"
jpayne@69: #include <vector>
jpayne@69: #include <algorithm>
jpayne@69: using namespace std;
jpayne@69: 
jpayne@69: //------------------------------------------------------------- Constants ----//
jpayne@69: 
jpayne@69: const char NUCMER_MISMATCH_CHAR = '^';
jpayne@69: const char NUCMER_MATCH_CHAR = ' ';
jpayne@69: const char PROMER_SIM_CHAR = '+';
jpayne@69: const char PROMER_MISMATCH_CHAR = ' ';
jpayne@69: 
jpayne@69: //-- Note: if coord exceeds LINE_PREFIX_LEN - 1 digits,
jpayne@69: //         increase these accordingly
jpayne@69: #define LINE_PREFIX_LEN 11
jpayne@69: #define PREFIX_FORMAT "%-10ld "
jpayne@69: 
jpayne@69: #define DEFAULT_SCREEN_WIDTH 60
jpayne@69: int Screen_Width = DEFAULT_SCREEN_WIDTH;
jpayne@69: 
jpayne@69: 
jpayne@69: 
jpayne@69: //------------------------------------------------------ Type Definitions ----//
jpayne@69: struct AlignStats
jpayne@69:      //-- Alignment statistics data structure
jpayne@69: {
jpayne@69:   long int sQ, eQ, sR, eR;              // start and end in Query and Reference
jpayne@69:                                         // relative to the directional strand
jpayne@69:   vector<long int> Delta;               // delta information
jpayne@69: };
jpayne@69: 
jpayne@69: 
jpayne@69: 
jpayne@69: struct sR_Sort
jpayne@69: //-- For sorting alignments by their sR coordinate
jpayne@69: {
jpayne@69:   bool operator( ) (const AlignStats & pA, const AlignStats & pB)
jpayne@69:   {
jpayne@69:     //-- sort sR
jpayne@69:     if ( pA.sR < pB.sR )
jpayne@69:       return true;
jpayne@69:     else
jpayne@69:       return false;
jpayne@69:   }
jpayne@69: };
jpayne@69: 
jpayne@69: 
jpayne@69: 
jpayne@69: struct sQ_Sort
jpayne@69: //-- For sorting alignments by their sQ coordinate
jpayne@69: {
jpayne@69:   bool operator( ) (const AlignStats & pA, const AlignStats & pB)
jpayne@69:   {
jpayne@69:     //-- sort sQ
jpayne@69:     if ( pA.sQ < pB.sQ )
jpayne@69:       return true;
jpayne@69:     else
jpayne@69:       return false;
jpayne@69:   }
jpayne@69: };
jpayne@69: 
jpayne@69: 
jpayne@69: 
jpayne@69: 
jpayne@69: //------------------------------------------------------ Global Variables ----//
jpayne@69: bool isSortByQuery = false;              // -q option
jpayne@69: bool isSortByReference = false;          // -r option
jpayne@69: 
jpayne@69: int DATA_TYPE = NUCMER_DATA;
jpayne@69: int MATRIX_TYPE = BLOSUM62;
jpayne@69: 
jpayne@69: char InputFileName [MAX_LINE];
jpayne@69: char RefFileName [MAX_LINE], QryFileName [MAX_LINE];
jpayne@69: 
jpayne@69: 
jpayne@69: 
jpayne@69: //------------------------------------------------- Function Declarations ----//
jpayne@69: long int toFwd
jpayne@69:      (long int coord, long int len, int frame);
jpayne@69: 
jpayne@69: void parseDelta
jpayne@69:      (vector<AlignStats> & Aligns, char * IdR, char * IdQ);
jpayne@69: 
jpayne@69: void printAlignments
jpayne@69:      (vector<AlignStats> Aligns, char * R, char * Q);
jpayne@69: 
jpayne@69: void printHelp
jpayne@69:      (const char * s);
jpayne@69: 
jpayne@69: void printUsage
jpayne@69:      (const char * s);
jpayne@69: 
jpayne@69: long int revC
jpayne@69:      (long int coord, long int len);
jpayne@69: 
jpayne@69: 
jpayne@69: 
jpayne@69: //-------------------------------------------------- Function Definitions ----//
jpayne@69: int main
jpayne@69:      (int argc, char ** argv)
jpayne@69: {
jpayne@69:   long int i;
jpayne@69: 
jpayne@69:   FILE * RefFile = NULL;
jpayne@69:   FILE * QryFile = NULL;
jpayne@69: 
jpayne@69:   vector<AlignStats> Aligns;
jpayne@69: 
jpayne@69:   char * R, * Q;
jpayne@69: 
jpayne@69:   long int InitSize = INIT_SIZE;
jpayne@69:   char Id [MAX_LINE], IdR [MAX_LINE], IdQ [MAX_LINE];
jpayne@69: 
jpayne@69:   //-- Parse the command line arguments
jpayne@69:   {
jpayne@69:     int ch, errflg = 0;
jpayne@69:     optarg = NULL;
jpayne@69: 
jpayne@69:     while ( !errflg  &&  ((ch = getopt
jpayne@69:                            (argc, argv, "hqrw:x:")) != EOF) )
jpayne@69:       switch (ch)
jpayne@69:         {
jpayne@69:         case 'h' :
jpayne@69: 	  printHelp (argv[0]);
jpayne@69: 	  exit (EXIT_SUCCESS);
jpayne@69:           break;
jpayne@69: 
jpayne@69: 	case 'q' :
jpayne@69: 	  isSortByQuery = true;
jpayne@69: 	  break;
jpayne@69: 
jpayne@69: 	case 'r' :
jpayne@69: 	  isSortByReference = true;
jpayne@69: 	  break;
jpayne@69: 
jpayne@69: 	case 'w' :
jpayne@69: 	  Screen_Width = atoi (optarg);
jpayne@69: 	  if ( Screen_Width <= LINE_PREFIX_LEN )
jpayne@69: 	    {
jpayne@69: 	      fprintf(stderr,
jpayne@69: 		      "WARNING: invalid screen width %d, using default\n",
jpayne@69: 		      DEFAULT_SCREEN_WIDTH);
jpayne@69: 	      Screen_Width = DEFAULT_SCREEN_WIDTH;
jpayne@69: 	    }
jpayne@69: 	  break;
jpayne@69: 
jpayne@69: 	case 'x' :
jpayne@69: 	  MATRIX_TYPE = atoi (optarg);
jpayne@69: 	  if ( MATRIX_TYPE < 1 || MATRIX_TYPE > 3 )
jpayne@69: 	    {
jpayne@69: 	      fprintf(stderr,
jpayne@69: 		      "WARNING: invalid matrix type %d, using default\n",
jpayne@69: 		      MATRIX_TYPE);
jpayne@69: 	      MATRIX_TYPE = BLOSUM62;
jpayne@69: 	    }
jpayne@69: 	  break;
jpayne@69: 
jpayne@69:         default :
jpayne@69:           errflg ++;
jpayne@69:         }
jpayne@69: 
jpayne@69:     if ( errflg > 0  ||  argc - optind != 3 )
jpayne@69:       {
jpayne@69:         printUsage (argv[0]);
jpayne@69:         exit (EXIT_FAILURE);
jpayne@69:       }
jpayne@69: 
jpayne@69:     if ( isSortByQuery  &&  isSortByReference )
jpayne@69:       fprintf (stderr,
jpayne@69:                "WARNING: both -r and -q were passed, -q ignored\n");
jpayne@69:   }
jpayne@69: 
jpayne@69:   strcpy (InputFileName, argv[optind ++]);
jpayne@69:   strcpy (IdR, argv[optind ++]);
jpayne@69:   strcpy (IdQ, argv[optind ++]);
jpayne@69: 
jpayne@69:   //-- Read in the alignment data
jpayne@69:   parseDelta (Aligns, IdR, IdQ);
jpayne@69: 
jpayne@69:   //-- Find, and read in the reference sequence
jpayne@69:   RefFile = File_Open (RefFileName, "r");
jpayne@69:   InitSize = INIT_SIZE;
jpayne@69:   R = (char *) Safe_malloc ( sizeof(char) * InitSize );
jpayne@69:   while ( Read_String (RefFile, R, InitSize, Id, FALSE) )
jpayne@69:     if ( strcmp (Id, IdR) == 0 )
jpayne@69:       break;
jpayne@69:   fclose (RefFile);
jpayne@69:   if ( strcmp (Id, IdR) != 0 )
jpayne@69:     {
jpayne@69:       fprintf(stderr,"ERROR: Could not find %s in the reference file\n", IdR);
jpayne@69:       exit (EXIT_FAILURE);
jpayne@69:     }
jpayne@69: 
jpayne@69: 
jpayne@69:   //-- Find, and read in the query sequence
jpayne@69:   QryFile = File_Open (QryFileName, "r");
jpayne@69:   InitSize = INIT_SIZE;
jpayne@69:   Q = (char *) Safe_malloc ( sizeof(char) * InitSize );
jpayne@69:   while ( Read_String (QryFile, Q, InitSize, Id, FALSE) )
jpayne@69:     if ( strcmp (Id, IdQ) == 0 )
jpayne@69:       break;
jpayne@69:   fclose (QryFile);
jpayne@69:   if ( strcmp (Id, IdQ) != 0 )
jpayne@69:     {
jpayne@69:       fprintf(stderr,"ERROR: Could not find %s in the query file\n", IdQ);
jpayne@69:       exit (EXIT_FAILURE);
jpayne@69:     }
jpayne@69: 
jpayne@69:   //-- Sort the alignment regions if user passed -r or -q option
jpayne@69:   if ( isSortByReference )
jpayne@69:     sort (Aligns.begin( ), Aligns.end( ), sR_Sort( ));
jpayne@69:   else if ( isSortByQuery )
jpayne@69:     sort (Aligns.begin( ), Aligns.end( ), sQ_Sort( ));
jpayne@69: 
jpayne@69: 
jpayne@69:   //-- Output the alignments to stdout
jpayne@69:   printf("%s %s\n\n", RefFileName, QryFileName);
jpayne@69:   for ( i = 0; i < Screen_Width; i ++ ) printf("=");
jpayne@69:   printf("\n-- Alignments between %s and %s\n\n", IdR, IdQ);
jpayne@69:   printAlignments (Aligns, R, Q);
jpayne@69:   printf("\n");
jpayne@69:   for ( i = 0; i < Screen_Width; i ++ ) printf("=");
jpayne@69:   printf("\n");
jpayne@69: 
jpayne@69:   return EXIT_SUCCESS;
jpayne@69: }
jpayne@69: 
jpayne@69: 
jpayne@69: 
jpayne@69: 
jpayne@69: long int toFwd
jpayne@69:      (long int coord, long int len, int frame)
jpayne@69: 
jpayne@69:      // Switch relative coordinate to reference forward DNA strand
jpayne@69: 
jpayne@69: {
jpayne@69:   long int newc = coord;
jpayne@69: 
jpayne@69:   if ( DATA_TYPE == PROMER_DATA )
jpayne@69:     newc = newc * 3 - (3 - labs(frame));
jpayne@69: 
jpayne@69:   if ( frame < 0 )
jpayne@69:     return revC ( newc, len );
jpayne@69:   else
jpayne@69:     return newc;
jpayne@69: }
jpayne@69: 
jpayne@69: 
jpayne@69: 
jpayne@69: 
jpayne@69: void parseDelta
jpayne@69:      (vector<AlignStats> & Aligns, char * IdR, char * IdQ)
jpayne@69: 
jpayne@69:      // Read in the alignments from the desired region
jpayne@69: 
jpayne@69: {
jpayne@69:   AlignStats aStats;                     //  single alignment region
jpayne@69:   bool found = false;
jpayne@69: 
jpayne@69:   DeltaReader_t dr;
jpayne@69:   dr.open (InputFileName);
jpayne@69:   DATA_TYPE = dr.getDataType( ) == NUCMER_STRING ?
jpayne@69:     NUCMER_DATA : PROMER_DATA;
jpayne@69:   strcpy (RefFileName, dr.getReferencePath( ).c_str( ));
jpayne@69:   strcpy (QryFileName, dr.getQueryPath( ).c_str( ));
jpayne@69: 
jpayne@69:   while ( dr.readNext( ) )
jpayne@69:     {
jpayne@69:       if ( dr.getRecord( ).idR == IdR  &&
jpayne@69: 	   dr.getRecord( ).idQ == IdQ )
jpayne@69: 	{
jpayne@69: 	  found = true;
jpayne@69: 	  break;
jpayne@69: 	}
jpayne@69:     }
jpayne@69:   if ( !found )
jpayne@69:     {
jpayne@69:       fprintf(stderr, "ERROR: Could not find any alignments for %s and %s\n",
jpayne@69: 	      IdR, IdQ);
jpayne@69:       exit (EXIT_FAILURE);
jpayne@69:     }
jpayne@69: 
jpayne@69:   for ( unsigned int i = 0; i < dr.getRecord( ).aligns.size( ); i ++ )
jpayne@69:     {
jpayne@69:       aStats.sR = dr.getRecord( ).aligns[i].sR;
jpayne@69:       aStats.eR = dr.getRecord( ).aligns[i].eR;
jpayne@69:       aStats.sQ = dr.getRecord( ).aligns[i].sQ;
jpayne@69:       aStats.eQ = dr.getRecord( ).aligns[i].eQ;
jpayne@69: 
jpayne@69:       aStats.Delta = dr.getRecord( ).aligns[i].deltas;
jpayne@69: 
jpayne@69:       //-- Add the new alignment
jpayne@69:       Aligns.push_back (aStats);
jpayne@69:     }
jpayne@69:   dr.close( );
jpayne@69: 
jpayne@69:   return;
jpayne@69: }
jpayne@69: 
jpayne@69: 
jpayne@69: 
jpayne@69: 
jpayne@69: void printAlignments
jpayne@69:      (vector<AlignStats> Aligns, char * R, char * Q)
jpayne@69: 
jpayne@69:      // Print the alignments to the screen
jpayne@69: 
jpayne@69: {
jpayne@69:   vector<AlignStats>::iterator Ap;
jpayne@69:   vector<long int>::iterator Dp;
jpayne@69: 
jpayne@69:   char * A[7] = {NULL, NULL, NULL, NULL, NULL, NULL, NULL};
jpayne@69:   char * B[7] = {NULL, NULL, NULL, NULL, NULL, NULL, NULL};
jpayne@69:   int Ai, Bi, i;
jpayne@69: 
jpayne@69:   char Buff1 [Screen_Width + 1],
jpayne@69:        Buff2 [Screen_Width + 1],
jpayne@69:        Buff3 [Screen_Width + 1];
jpayne@69: 
jpayne@69:   int Sign;
jpayne@69:   long int Delta;
jpayne@69:   long int Total, Errors, Remain;
jpayne@69:   long int Pos;
jpayne@69: 
jpayne@69:   long int sR, eR, sQ, eQ;
jpayne@69:   long int Apos, Bpos;
jpayne@69:   long int SeqLenR, SeqLenQ;
jpayne@69:   int frameR, frameQ;
jpayne@69: 
jpayne@69:   for ( i = 0; i < LINE_PREFIX_LEN; i ++ )
jpayne@69:     Buff3[i] = ' ';
jpayne@69: 
jpayne@69:   SeqLenR = strlen (R + 1);
jpayne@69:   SeqLenQ = strlen (Q + 1);
jpayne@69: 
jpayne@69:   if ( DATA_TYPE == NUCMER_DATA )
jpayne@69:     {
jpayne@69:       A[1] = R;
jpayne@69:       A[4] = (char *) Safe_malloc ( sizeof(char) * (SeqLenR + 2) );
jpayne@69:       strcpy ( A[4] + 1, A[1] + 1 );
jpayne@69:       A[4][0] = '\0';
jpayne@69:       Reverse_Complement ( A[4], 1, SeqLenR );
jpayne@69: 
jpayne@69:       B[1] = Q;
jpayne@69:       B[4] = (char *) Safe_malloc ( sizeof(char) * (SeqLenQ + 2) );
jpayne@69:       strcpy ( B[4] + 1, B[1] + 1 );
jpayne@69:       B[4][0] = '\0';
jpayne@69:       Reverse_Complement ( B[4], 1, SeqLenQ );
jpayne@69:     }
jpayne@69: 
jpayne@69:   for ( Ap = Aligns.begin( ); Ap < Aligns.end( ); Ap ++ )
jpayne@69:     {
jpayne@69:       sR = Ap->sR;
jpayne@69:       eR = Ap->eR;
jpayne@69:       sQ = Ap->sQ;
jpayne@69:       eQ = Ap->eQ;
jpayne@69: 
jpayne@69:       //-- Get the coords and frame right
jpayne@69:       frameR = 1;
jpayne@69:       if ( sR > eR )
jpayne@69: 	{
jpayne@69: 	  sR = revC (sR, SeqLenR);
jpayne@69: 	  eR = revC (eR, SeqLenR);
jpayne@69: 	  frameR += 3;
jpayne@69: 	}
jpayne@69:       frameQ = 1;
jpayne@69:       if ( sQ > eQ )
jpayne@69: 	{
jpayne@69: 	  sQ = revC (sQ, SeqLenQ);
jpayne@69: 	  eQ = revC (eQ, SeqLenQ);
jpayne@69: 	  frameQ += 3;
jpayne@69: 	}
jpayne@69: 
jpayne@69:       if ( DATA_TYPE == PROMER_DATA )
jpayne@69: 	{
jpayne@69: 	  frameR += (sR + 2) % 3;
jpayne@69: 	  frameQ += (sQ + 2) % 3;
jpayne@69: 
jpayne@69: 	  //-- Translate the coordinates from DNA to Amino Acid
jpayne@69: 	  //   remeber that eR and eQ point to the last nucleotide in the codon
jpayne@69: 	  sR = (sR + 2) / 3;
jpayne@69: 	  eR = eR / 3;
jpayne@69: 	  sQ = (sQ + 2) / 3;
jpayne@69: 	  eQ = eQ / 3;
jpayne@69: 	}
jpayne@69:       Ai = frameR;
jpayne@69:       Bi = frameQ;
jpayne@69:       if ( frameR > 3 )
jpayne@69: 	frameR = -(frameR - 3);
jpayne@69:       if ( frameQ > 3 )
jpayne@69: 	frameQ = -(frameQ - 3);
jpayne@69: 
jpayne@69:       //-- Get the sequence
jpayne@69:       if ( A[Ai] == NULL )
jpayne@69: 	{
jpayne@69: 	  assert ( DATA_TYPE == PROMER_DATA );
jpayne@69: 	  A[Ai] = (char *) Safe_malloc ( sizeof(char) * ( SeqLenR / 3 + 2 ) );
jpayne@69: 	  A[Ai][0] = '\0';
jpayne@69: 	  Translate_DNA ( R, A[Ai], Ai );
jpayne@69: 	}
jpayne@69:       if ( B[Bi] == NULL )
jpayne@69: 	{
jpayne@69: 	  assert ( DATA_TYPE == PROMER_DATA );
jpayne@69: 	  B[Bi] = (char *) Safe_malloc ( sizeof(char) * ( SeqLenQ / 3 + 2 ) );
jpayne@69: 	  B[Bi][0] = '\0';
jpayne@69: 	  Translate_DNA ( Q, B[Bi], Bi );
jpayne@69: 	}
jpayne@69: 
jpayne@69: 
jpayne@69:       //-- Generate the alignment
jpayne@69:       printf("-- BEGIN alignment [ %s%d %ld - %ld | %s%d %ld - %ld ]\n\n\n",
jpayne@69: 	     frameR > 0 ? "+" : "-", abs(frameR), Ap->sR, Ap->eR,
jpayne@69: 	     frameQ > 0 ? "+" : "-", abs(frameQ), Ap->sQ, Ap->eQ);
jpayne@69: 
jpayne@69:       Apos = sR;
jpayne@69:       Bpos = sQ;
jpayne@69: 
jpayne@69:       Errors = 0;
jpayne@69:       Total = 0;
jpayne@69:       Remain = eR - sR + 1;
jpayne@69: 
jpayne@69:       sprintf(Buff1, PREFIX_FORMAT, toFwd (Apos, SeqLenR, frameR));
jpayne@69:       sprintf(Buff2, PREFIX_FORMAT, toFwd (Bpos, SeqLenQ, frameQ));
jpayne@69:       Pos = LINE_PREFIX_LEN;
jpayne@69: 
jpayne@69:       for ( Dp = Ap->Delta.begin( );
jpayne@69: 	    Dp < Ap->Delta.end( ) &&
jpayne@69: 	    *Dp != 0; Dp ++ )
jpayne@69: 	{
jpayne@69: 	  Delta = *Dp;
jpayne@69: 	  Sign = Delta > 0 ? 1 : -1;
jpayne@69: 	  Delta = labs ( Delta );
jpayne@69: 
jpayne@69: 
jpayne@69: 	  //-- For all the bases before the next indel
jpayne@69: 	  for ( i = 1; i < Delta; i ++ )
jpayne@69: 	    {
jpayne@69: 	      if ( Pos >= Screen_Width )
jpayne@69: 		{
jpayne@69: 		  Buff1[Pos] = Buff2[Pos] = Buff3[Pos] = '\0';
jpayne@69: 		  if ( DATA_TYPE == NUCMER_DATA )
jpayne@69: 		    printf("%s\n%s\n%s\n\n", Buff1, Buff2, Buff3);
jpayne@69: 		  else
jpayne@69: 		    printf("%s\n%s\n%s\n\n", Buff1, Buff3, Buff2);
jpayne@69: 		  sprintf(Buff1, PREFIX_FORMAT, toFwd (Apos, SeqLenR, frameR));
jpayne@69: 		  sprintf(Buff2, PREFIX_FORMAT, toFwd (Bpos, SeqLenQ, frameQ));
jpayne@69: 		  Pos = LINE_PREFIX_LEN;
jpayne@69: 		}
jpayne@69: 
jpayne@69: 	      if ( DATA_TYPE == NUCMER_DATA )
jpayne@69: 		Buff3[Pos] = A[Ai][Apos] == B[Bi][Bpos] ?
jpayne@69: 		  NUCMER_MATCH_CHAR : NUCMER_MISMATCH_CHAR;
jpayne@69: 	      else if ( A[Ai][Apos] == B[Bi][Bpos] )
jpayne@69: 		Buff3[Pos] = A[Ai][Apos];
jpayne@69: 	      else
jpayne@69: 		Buff3[Pos] = MATCH_SCORE
jpayne@69: 		  [MATRIX_TYPE]
jpayne@69: 		  [toupper(A[Ai][Apos]) - 'A']
jpayne@69: 		  [toupper(B[Bi][Bpos]) - 'A'] > 0 ?
jpayne@69: 		  PROMER_SIM_CHAR : PROMER_MISMATCH_CHAR;
jpayne@69: 	      Buff1[Pos] = A[Ai][Apos ++];
jpayne@69: 	      Buff2[Pos ++] = B[Bi][Bpos ++];
jpayne@69: 	    }
jpayne@69: 
jpayne@69: 
jpayne@69: 	  //-- For the indel
jpayne@69: 	  Remain -= i - 1;
jpayne@69: 
jpayne@69: 	  if ( Pos >= Screen_Width )
jpayne@69: 	    {
jpayne@69: 	      Buff1[Pos] = Buff2[Pos] = Buff3[Pos] = '\0';
jpayne@69: 	      if ( DATA_TYPE == NUCMER_DATA )
jpayne@69: 		printf("%s\n%s\n%s\n\n", Buff1, Buff2, Buff3);
jpayne@69: 	      else
jpayne@69: 		printf("%s\n%s\n%s\n\n", Buff1, Buff3, Buff2);
jpayne@69: 	      sprintf(Buff1, PREFIX_FORMAT, toFwd (Apos, SeqLenR, frameR));
jpayne@69: 	      sprintf(Buff2, PREFIX_FORMAT, toFwd (Bpos, SeqLenQ, frameQ));
jpayne@69: 	      Pos = LINE_PREFIX_LEN;
jpayne@69: 	    }
jpayne@69: 
jpayne@69: 	  if ( Sign == 1 )
jpayne@69: 	    {
jpayne@69: 	      if ( DATA_TYPE == NUCMER_DATA )
jpayne@69: 		Buff3[Pos] = NUCMER_MISMATCH_CHAR;
jpayne@69: 	      else
jpayne@69: 		Buff3[Pos] = PROMER_MISMATCH_CHAR;
jpayne@69: 	      Buff1[Pos] = A[Ai][Apos ++];
jpayne@69: 	      Buff2[Pos ++] = '.';
jpayne@69: 	      Remain --;
jpayne@69: 	    }
jpayne@69: 	  else
jpayne@69: 	    {
jpayne@69: 	      if ( DATA_TYPE == NUCMER_DATA )
jpayne@69: 		Buff3[Pos] = NUCMER_MISMATCH_CHAR;
jpayne@69: 	      else
jpayne@69: 		Buff3[Pos] = PROMER_MISMATCH_CHAR;
jpayne@69: 	      Buff1[Pos] = '.';
jpayne@69: 	      Buff2[Pos ++] = B[Bi][Bpos ++];
jpayne@69: 	      Total ++;
jpayne@69: 	    }
jpayne@69: 	}
jpayne@69: 
jpayne@69: 
jpayne@69:       //-- For all the bases remaining after the last indel
jpayne@69:       for ( i = 0; i < Remain; i ++ )
jpayne@69: 	{
jpayne@69: 	  if ( Pos >= Screen_Width )
jpayne@69: 	    {
jpayne@69: 	      Buff1[Pos] = Buff2[Pos] = Buff3[Pos] = '\0';
jpayne@69: 	      if ( DATA_TYPE == NUCMER_DATA )
jpayne@69: 		printf("%s\n%s\n%s\n\n", Buff1, Buff2, Buff3);
jpayne@69: 	      else
jpayne@69: 		printf("%s\n%s\n%s\n\n", Buff1, Buff3, Buff2);
jpayne@69: 	      sprintf(Buff1, PREFIX_FORMAT, toFwd (Apos, SeqLenR, frameR));
jpayne@69: 	      sprintf(Buff2, PREFIX_FORMAT, toFwd (Bpos, SeqLenQ, frameQ));
jpayne@69: 	      Pos = LINE_PREFIX_LEN;
jpayne@69: 	    }
jpayne@69: 	  
jpayne@69: 	  if ( DATA_TYPE == NUCMER_DATA )
jpayne@69: 	    Buff3[Pos] = A[Ai][Apos] == B[Bi][Bpos] ?
jpayne@69: 	      NUCMER_MATCH_CHAR : NUCMER_MISMATCH_CHAR;
jpayne@69: 	  else if ( A[Ai][Apos] == B[Bi][Bpos] )
jpayne@69: 	    Buff3[Pos] = A[Ai][Apos];
jpayne@69: 	  else
jpayne@69: 	    Buff3[Pos] = MATCH_SCORE
jpayne@69: 	      [MATRIX_TYPE]
jpayne@69: 	      [toupper(A[Ai][Apos]) - 'A']
jpayne@69: 	      [toupper(B[Bi][Bpos]) - 'A'] > 0 ?
jpayne@69: 	      PROMER_SIM_CHAR : PROMER_MISMATCH_CHAR;
jpayne@69: 	  Buff1[Pos] = A[Ai][Apos ++];
jpayne@69: 	  Buff2[Pos ++] = B[Bi][Bpos ++];
jpayne@69: 	}
jpayne@69: 
jpayne@69: 
jpayne@69:       //-- For the remaining buffered output
jpayne@69:       if ( Pos > LINE_PREFIX_LEN )
jpayne@69: 	{
jpayne@69: 	  Buff1[Pos] = Buff2[Pos] = Buff3[Pos] = '\0';
jpayne@69: 	  if ( DATA_TYPE == NUCMER_DATA )
jpayne@69: 	    printf("%s\n%s\n%s\n\n", Buff1, Buff2, Buff3);
jpayne@69: 	  else
jpayne@69: 	    printf("%s\n%s\n%s\n\n", Buff1, Buff3, Buff2);
jpayne@69: 	  sprintf(Buff1, PREFIX_FORMAT, toFwd (Apos, SeqLenR, frameR));
jpayne@69: 	  sprintf(Buff2, PREFIX_FORMAT, toFwd (Bpos, SeqLenQ, frameQ));
jpayne@69: 	  Pos = LINE_PREFIX_LEN;
jpayne@69: 	}
jpayne@69: 
jpayne@69:       printf("\n--   END alignment [ %s%d %ld - %ld | %s%d %ld - %ld ]\n",
jpayne@69: 	     frameR > 0 ? "+" : "-", abs(frameR), Ap->sR, Ap->eR,
jpayne@69: 	     frameQ > 0 ? "+" : "-", abs(frameQ), Ap->sQ, Ap->eQ);
jpayne@69:     }
jpayne@69: 
jpayne@69:   //-- Free the sequences, except for the originals
jpayne@69:   for ( i = 0; i < 7; i ++ )
jpayne@69:     {
jpayne@69:       if ( (DATA_TYPE != NUCMER_DATA || i != 1)  &&  A[i] != NULL )
jpayne@69: 	free ( A[i] );
jpayne@69:       if ( (DATA_TYPE != NUCMER_DATA || i != 1)  &&  B[i] != NULL )
jpayne@69: 	free ( B[i] );
jpayne@69:     }
jpayne@69: 
jpayne@69:   return;
jpayne@69: }
jpayne@69: 
jpayne@69: 
jpayne@69: 
jpayne@69: 
jpayne@69: void printHelp
jpayne@69:      (const char * s)
jpayne@69: 
jpayne@69:       //  Display the program's help information to stderr
jpayne@69: 
jpayne@69: {
jpayne@69:   fprintf (stderr,
jpayne@69:            "\nUSAGE: %s  [options]  <deltafile>  <ref ID>  <qry ID>\n\n", s);
jpayne@69:   fprintf (stderr,
jpayne@69:        "-h            Display help information\n"
jpayne@69:        "-q            Sort alignments by the query start coordinate\n"
jpayne@69:        "-r            Sort alignments by the reference start coordinate\n"
jpayne@69:        "-w int        Set the screen width - default is 60\n"
jpayne@69:        "-x int        Set the matrix type - default is 2 (BLOSUM 62),\n"
jpayne@69:        "              other options include 1 (BLOSUM 45) and 3 (BLOSUM 80)\n"
jpayne@69:        "              note: only has effect on amino acid alignments\n\n");
jpayne@69:   fprintf (stderr,
jpayne@69:        "  Input is the .delta output of either the \"nucmer\" or the\n"
jpayne@69:        "\"promer\" program passed on the command line.\n"
jpayne@69:        "  Output is to stdout, and consists of all the alignments between the\n"
jpayne@69:        "query and reference sequences identified on the command line.\n"
jpayne@69:        "  NOTE: No sorting is done by default, therefore the alignments\n"
jpayne@69:        "will be ordered as found in the <deltafile> input.\n\n");
jpayne@69:   return;
jpayne@69: }
jpayne@69: 
jpayne@69: 
jpayne@69: 
jpayne@69: 
jpayne@69: void printUsage
jpayne@69:      (const char * s)
jpayne@69: 
jpayne@69:       //  Display the program's usage information to stderr.
jpayne@69: 
jpayne@69: {
jpayne@69:   fprintf (stderr,
jpayne@69:            "\nUSAGE: %s  [options]  <deltafile>  <ref ID>  <qry ID>\n\n", s);
jpayne@69:   fprintf (stderr, "Try '%s -h' for more information.\n", s);
jpayne@69:   return;
jpayne@69: }
jpayne@69: 
jpayne@69: 
jpayne@69: 
jpayne@69: 
jpayne@69: long int revC
jpayne@69:      (long int coord, long int len)
jpayne@69: {
jpayne@69:   return len - coord + 1;
jpayne@69: }