/*=====================================================================*/
/*(C) Copyright 1991-9 by Fred Hutchinson Cancer Research Center         */
/*        multimat.c  reads several BLIMPS MATRIX search output files
	   and merges the results.  Writes stats to multimat.dat.
	    eg:  multimat 10 blocks.dat PS00094.lis PS00094*.hom or
		 multimat 15 none none PS00094A.hom PS00094B.hom ... (up to 26)

	   1st arg = number of hits to report
	   2nd arg = name of blocks database (use non-existent file name
		     to skip it, if file is found, will open it & compare
		     search results with blocks of same name)
           3rd arg = name of list of true positive sequences (use
                     non-existent file name to omit it) If file is found
                     will only report hits NOT in file. 
	   4rd arg = name of matrix search output files, etc.
			EG PS00094A, PS00094B, PS00094C

	Name restrictions:
		AC name in the blocks in the blocks database (2nd arg)
			must be unique and <= 8 characters long.
			The blocks database must be sorted in AC order.
			It is assumed that all of the blocks from the
			same family have all the same AC except for the
			last character.
		Search result file names (4th arg) must match the AC names of
			the blocks that were the queries for the search.
		Sequence names in the blocks in the blocks db (2nd arg)
			and in the list of true positives (3rd arg)
			must be unique, and must match the sequence
			names in the database searched. It is assumed that
			searches are for blocks belonging to the same
			family.
--------------------------------------------------------------------------
     8/1/91  J. Henikoff

   5/15/99 1. Format changes for Blimps 3.2.6
   6/26/99 1. Sequence name length from 10 to 18
   6/11/00 1. AC length from 7 to 9
  12/23/06 1. Sequence name length from 18 to 20
====================================================================*/
#include "motifj.h"

#define SWISS 69113		/* Swiss35 */
#define NSCORE 5000		/* Maximum # of scores in .hom file */
#define MAXHOM 26		/* Maximum # of .hom files */
#define MAXHIT 10		/* Maximum # of hits to report */
#define MAXMEM 64000		/* Maximum bytes for DOS array */
#define SEEK_SET 0
#define PROTEIN 1
#define DNA 3
#define MAXMAP 60		/* Maximum map width in characters */
#define MAX_WIDTH 100		/* Maximum block width expected */

struct hom {			/* search results structure */
   char ac[MAXAC+1];		/* block name, eg. PS00094A */
   int norm;			/* normalized score */
   int max_norm;		/* max normalized score for this group */
   int min_rank;		/* minimum rank for this group of blocks */
   int rank;			/* rank of this result in search */
   int frame;			/* frame of alignment */
   int strength;		/* block strengh */
   int score;			/* patmat score */
   int offset;			/* offset of alignment */
   double seqlen;			/* db sequence length */
   int width;			/* width of block */
   char title[30];		/* block description */
   char seq_id[25];		/* patmat sequence id */
   char aa[MAX_WIDTH]; 	/* alignment to block */
   int map_flag;		/* used by map_blocks */
   int tp;			/* true positive flag, 1 if tp, 2 if in block*/
};

struct best_block {			/* block structure */
   char ac[MAXAC+1];         		/* accession number */
   int nseq, width, strength;		/* #seqs, width, strength */
   int s995;				/* 99.5 percentile score */
   int minprev, maxprev;		/* distances from prev block */
   char name[MAXSEQS][SNAMELEN+5];		/* name of seq */
   int offset[MAXSEQS];			/* offset of seq */
   char aa[MAXSEQS][MAX_WIDTH];	/* aas for seq */
   int cluster[MAXSEQS];		/* cluster # for seq */
   int ncluster[MAXSEQS];		/* #seqs in same cluster */
   long dat_pos;			/* position of block in blocks.dat*/
   int query;				/* score if last seq is a query */
   int rank;				/* rank if last seq is a query */
   int anchor;				/* anchor block flag */
   struct best_block *next_block;	/* next block in path */
};

/*---- Functions also in blksort.c -----*/
void print_blurb();
void check_dat();
struct best_block *read_block();
void fill_block();
int add_query();
double hypergeo();
int distance_okay();
int prev_dist();
int distance();
void map_blocks();
void align_blocks();
int closest_seq();
int tempcmp();
int read_hom();
int check_tp();
void show_hom();
int id_cmp();
int rank_cmp();
int idnorm_cmp();
int norm_cmp();
struct best_block *get_blocks();
void consensus();
int compute_loc();
void re_normalize();
/*--------------------Routines from motmisc.o------------------------*/
struct db_id *makedbid();
int get_ids();
struct split_name *split_names();
void kr_itoa();

/*------------Global variables -------------*/
char Version[12] = " 6/11/00.1";/* Version number */
int Ask;			/* Interactive flag */
char AC[30];				/* Global variables */
int ACLen;				/* Length of AC */
int NScore;				/* Max# scores for DOS */
int LisSeq;				/* #seqs in TP list file */
int FragSeq;
int BlkSeq;
int SeqType;			/* DNA or PROTEIN database searched */
char Query[SNAMELEN];		/* Query file name */
double DBLen;			/* DB sequence length */
int NBlock;			/* # blocks in family */
double NSeq;			/* # sequences searched */
double NRank;			/* # of ranks assigned in search */
int DBType;			/* =1 if amino acid, =3 if nucleotide */
int MaxHit;			/* # hits to report */
FILE *Fout;			/* multimat.lis statistics file (not TPs) */
FILE *Fdat;			/* multimat.dat statistics file (TPs) */
FILE *Fmis;			/* multimat.mis file (missed TPs)     */
FILE *Ffnd;			/* multimat.fnd file (missed TPs)     */
/*======================================================================*/
int main(argc, argv)
int argc;
char *argv[];
{
   FILE *fhom, *fblk, *flis;
   char homfile[MAXHOM][FNAMELEN], *ptr, *ptr1, ctemp[FNAMELEN];
   char datfile[FNAMELEN], lisfile[FNAMELEN];
   int totscores, i, nhom;
   struct hom *temp;
   struct db_id *ids, *did;
   struct best_block *blocks;

   printf("MULTIMAT: (C) Copyright 1991 by Fred Hutchinson Cancer ");
   printf("Research Center\n");
   printf("Version %s\n", Version);
   if (argc <= 4)
   {
      printf("USAGE:   multimat <nhits> <blocks> <list> <blimps>\n");
      printf("            <nhits>  = Maximum # of hits to report\n");
      printf("            <blocks> = File containing blocks for comparison\n");
      printf("            <list>   = File containing list of TP sequences\n");
      printf("            <blimps> = Blimps search result files, ");
      printf(" blocks vs sequence db\n");
   }
/*   Fout = fopen("multimat.lis", "a");	 */
   Fdat = fopen("multimat.dat", "a");		/* statistics file */

/*--------------arg 1: Number of hits to report---------------------*/
   MaxHit = MAXHIT;
   ctemp[0] = '\0';
   if (argc > 1) strcpy(ctemp, argv[1]);
   else
   {
      printf("\nEnter maximum number of hits to report [%d]: ", MaxHit);
      gets(ctemp);
   }
   if (strlen(ctemp)) MaxHit = atoi(ctemp);
   if (MaxHit < 1 || MaxHit > NSCORE) MaxHit = MAXHIT;
/*------------- arg 2:  blocks database file -----------------------------*/
   if (argc > 2)
      strcpy(datfile, argv[2]);
   else
   {
      printf("\nEnter name of blocks database:\n");
      gets(datfile);
   }
   if ( (fblk=fopen(datfile, "r")) == NULL)
   {
      printf("\nCannot open file %s", datfile);
      printf(";\nsearch results will not be compared with blocks database.");
   }
   else
   {
      printf("\nSearch results will be compared with %s", datfile);
   }

/*------------- arg 3:  .lis file ----------------------------------------*/
   if (argc > 3)
      strcpy(lisfile, argv[3]);
   else
   {
      printf("\nEnter name of file containing list of true positives:\n");
      gets(lisfile);
   }
   if ( (flis=fopen(lisfile, "r")) == NULL)
   {
      printf("\nCannot open file %s;", lisfile);
    printf("\nsearch results will not be compared with true positive list.\n");
      LisSeq = FragSeq = BlkSeq = 0;  ids = NULL;
   }
   else
   {
      printf("\nSearch results will be compared with %s;", lisfile);
      printf("\nonly hits not in this list will be reported\n");
      ids = makedbid();
      printf("%d IDs in %s\n", (LisSeq=get_ids(flis, ids)), lisfile);
      fclose(flis); 
      /*----------- Count the true positive sequences -------------------*/
      LisSeq = FragSeq = BlkSeq = 0;
      did = ids->next;
      while (did != NULL) 
      {
         did->found = NO;
         if (!did->block && !did->frag)  	LisSeq++;
         if (did->block)           BlkSeq++;
         if (did->frag)            FragSeq++;
         did = did->next;
      }
      Fmis = fopen("multimat.mis", "w");		/* statistics file */
      Ffnd = fopen("multimat.fnd", "w");		/* statistics file */
   }

/*------------- args 4-:  .hom files -------------------------------*/
   nhom = 0;
   if (argc > MAXHOM+4) argc=MAXHOM+4;
   if (argc > 4)
   {
      for (i=4; i<argc; i++)
         strcpy(homfile[nhom++], argv[i]);
   }
   else
   {
      do
      {
	 printf("\nEnter name of BLIMPS results file: ");
	 gets(homfile[nhom++]);
      } while (nhom <= MAXHOM && strlen(homfile[nhom-1]));
      nhom--;
   }
   if (nhom < 1)
   {
      printf("\nNo results files specified.\n");
      exit(-1);
   }

/*----------Create an array big enough for NSCORE scores per file------*/
   NScore=NSCORE;
/*DOS  For DOS, have to restrict max size of array to 64k
   ltemp = MAXMEM; ltemp = ltemp/(nhom*sizeof(struct hom));
   NScore = (int) ltemp;
DOS*/
   temp = (struct hom *) malloc(nhom * NScore * sizeof(struct hom));
   if (temp==NULL)
   {
      printf("\n\nNOT ENOUGH MEMORY!\n");
      exit(-1);
   }

/*-------------Now read the files into the array for sorting-----------*/
   totscores = 0;
   DBType = PROTEIN;     /* BLIMPS always reports translated length */
   for (i=0; i<nhom; i++)
   {
       if ( (fhom=fopen(homfile[i], "r")) == NULL)
       {
	    printf("\nCannot open file %s\n", homfile[i]);
	    exit(-1);
       }
       else printf("\nReading %s...", homfile[i]);
       ptr = strrchr(homfile[i],'/');		/* look for file name */
       if (ptr != NULL)  ptr1 = strtok(ptr+1, ".");
       else  ptr1 = strtok(homfile[i], ".");
       strcpy(AC, ptr1);
       NScore=read_hom(fhom, temp, totscores, ids);
       printf("\n%d scores read from %s\n", NScore, homfile[i]);
       totscores += NScore;
       fclose(fhom);
   }
/*   AC is assumed to be a family name; all but the last character of
     the file name, maximum of 9 characters  */
   if (nhom > 1) AC[strlen(AC)-1] = '\0';
   if ((int) strlen(AC) > 9)
   {
      printf("\nWARNING: AC family name truncated to 9 characters. ");
      AC[9] = '\0';
   }
   printf("AC family name = %s\n", AC);
   ACLen = strlen(AC);

/*------------------Make a list of the blocks in the database----*/
   /*    For this AC => all but last character of AC must be the same for */
   /*    all block queries */
   if (fblk != NULL)
   {
      blocks = get_blocks(fblk, AC);
      re_normalize(temp, totscores, blocks);
      fclose(fblk);
   }
   else blocks = NULL;

/*-----------------Sort and present the results -----------------------*/
   show_hom(temp, totscores, blocks);

/*------------Write a file of the tps not found ------------------------*/
   if (ids != NULL)
   {
      did = ids->next;
      while (did != NULL) 
      {
         if (!did->found)
         {
            fprintf(Fmis, "%s", did->entry);
            if (did->block) fprintf(Fmis, "\tBLOCK");
            if (did->frag) fprintf(Fmis, "\tFRAGMENT");
            fprintf(Fmis, "\n");
         }
         else
         {
            fprintf(Ffnd, "%s", did->entry);
            if (did->block) fprintf(Ffnd, "\tBLOCK");
            if (did->frag) fprintf(Ffnd, "\tFRAGMENT");
            fprintf(Ffnd, "\n");
         }
         did = did->next;
      }
   }
/*---------------------------------------------------------------------*/
   if (Fout != NULL) fclose(Fout);
   if (Fdat != NULL) fclose(Fdat);
   if (Fmis != NULL) fclose(Fmis);
   if (Ffnd != NULL) fclose(Ffnd);
   printf("\n");
   exit(0);
}  /*  end of main */
/*============================================================================
BLIMPS (BLocks IMProved Searcher)   Version 3.2.5  1998/05
(C) Copyright 1993-2000, Fred Hutchinson Cancer Research Center



Block File: /howard/jorja/seqs/PACA.blk

Target File (s) : stuff.dna

Records Searched:   18

Scores Done:        108

Alignments Done:    32866

AC#                      Description                                                Score RF  AA# Length
pa01l_141.s1         0   CHROMAT_FILE: pa01l_141.s1 PHD_FILE: pa01l_141.s1.phd.1 TI 2081 -1   116    260 rYArlvKemsEkvqfiyithnkiaMEMAdqlmgvTmhEpgcsrlVavDveeav

NOTE: Uses hard-coded positions on Blimps output line, mainly because there
      are spaces in the Description.
=============================================================================*/
int read_hom(fhom, temp, first, tps)
FILE *fhom;
struct hom *temp;
int first;
struct db_id *tps;
{
   int t, i, v324, v326, minlen, offscore, offframe, offoffset, offlen;
   int offtitle;
   char line[MAXLINE], ctemp[30], *ptr;

   SeqType = PROTEIN;		/* default database type */
   t = first;    /* first available position in temp array */
   v324 = NO;	/* format changed with Blimps 3.2.4  */
   v326 = YES;  /* format changed again with Blimps 3.2.6 */
   minlen = 105; offscore = 84; offframe = 89; offoffset = 92; offlen = 98;
   offtitle = 22;
   while (!feof(fhom) && fgets(line, MAXLINE, fhom) != NULL)
   {
      if (strstr(line, "Version 3.2.3") ||
          strstr(line, "Version 3.2.2") || strstr(line, "Version 3.2.1")   )
      {  
         v324 = v326 = NO;
         minlen = 96; offscore = 75; offframe = 80; offoffset = 83; offlen = 89;
         offtitle = 14;
      }
      else if ( strstr(line, "Version 3.2.4") || strstr(line, "Version 3.2.5") )
      {  
         v324 = YES; v326 = NO;
         minlen =105; offscore = 84; offframe = 89; offoffset = 92; offlen = 98;
         offtitle = 24;
      }
      else if (strstr(line, "Records Searched:"))
      {
         ptr = strtok(line, ":"); ptr = strtok(NULL, "\t\r\n");
	 NSeq = atof(ptr);
      }
      else if (strstr(line, "Scores Done:"))
      {
         ptr = strtok(line, ":"); ptr = strtok(NULL, "\t\r\n");
	 NRank = atof(ptr);
      }
      else if ((int) strlen(line) > minlen && (t-first) < NScore)
      {
	 strcpy(temp[t].ac, AC);
         if (v324 || v326)
         {
	    strncpy(temp[t].seq_id, &line[0], 20); temp[t].seq_id[20]='\0';
         }
         else
         {
	    strncpy(temp[t].seq_id, &line[0], 12); temp[t].seq_id[12]='\0';
         }
	 strncpy(temp[t].title, &line[offtitle], 19); 
         temp[t].title[19] = '\0';
	 strncpy(ctemp, &line[offscore], 4); ctemp[4]='\0';
	 temp[t].score = atoi(ctemp);
	 strncpy(ctemp, &line[offframe], 2); ctemp[2]='\0';
	 temp[t].frame = atoi(ctemp);
	 if (temp[t].frame != 0) SeqType = DNA;
	 /*   BLIMPS frames are only non-zero for DNA databases */
	 if (temp[t].frame != 0) DBType = DNA;
	 strncpy(ctemp, &line[offoffset], 5); ctemp[5]='\0';
	 temp[t].offset = atoi(ctemp);
	 strncpy(ctemp, &line[offlen], 6); ctemp[6]='\0';
	 temp[t].seqlen = atof(ctemp);
	 temp[t].width = 0; i=minlen;

	 /*  truncates wide blocks */
	 while (line[i] != '\n' && temp[t].width < MAX_WIDTH) 
         {
            temp[t].aa[temp[t].width] = line[i++];
            temp[t].width += 1;
         }
         temp[t].aa[temp[t].width] = '\0';
         strcpy(ctemp, temp[t].seq_id);
         if (tps != NULL)
            temp[t].tp = check_tp(ctemp, tps);
	 temp[t].min_rank = temp[t].rank = t-first;
	 temp[t].map_flag = NO;
	 t++;
      }
   }
   /*    Compute normalized score */
   for (i = first; i < t; i++)
   {
      temp[i].norm = round((float) 1000. * temp[i].score / temp[t-1].score);
      temp[i].max_norm = temp[i].norm;
   }
   return(t-first);
}  /* end of read_hom */
/*=====================================================================
    See if sequence id is in the list of true positives
	Write info to multimat.dat if it's a TP not in the blocks
========================================================================*/
int check_tp(id, tps)
char *id;
struct db_id *tps;
{
   struct db_id *did;
   
   did = tps->next;
   while (did != NULL)
   {
      if (strncmp(id, did->entry, strlen(did->entry)) == 0)
      {
	 did->found = YES;
         if (did->block) return(2);
         else return(YES);
      }
      did = did->next;
   }
   return(NO);
}  /* end of check_tp */
/*=====================================================================*/
void show_hom(temp, tot, blocks)
struct hom *temp;
int tot;
struct best_block *blocks;
{
   int i, t, nhit, ntp, itn, itp, nblock, save_t, pearson, tp999, tn999;
   int tn[NSCORE], tp[NSCORE], rank, nfrag;
   int save_norm, save_strand, strand;
   long locfirst, loclast;
   double roc;
   char save_id[25];

   /*   Sort by strand, seq id, normalized score (d)  */
   qsort(temp, tot, sizeof(struct hom), idnorm_cmp);
   /*   Propagate maximum normalized score to all blocks for a group */
   strcpy(save_id, temp[0].seq_id); save_norm = temp[0].max_norm;
   if (temp[0].frame < 0) save_strand = -1;
   else                   save_strand = 1;
   save_t = 0;
   for (t=1; t<tot; t++)
   {
      if (temp[t].frame < 0) strand = -1; else strand = 1;
      if (save_strand == strand &&
          strcmp(temp[t].seq_id, save_id)==0)
      {
	 temp[t].max_norm = save_norm;
      }
      else
      {

	 strcpy(save_id, temp[t].seq_id);
         save_norm = temp[t].norm;
	 save_t = t;
         save_strand = strand;
      }
   }
   /*   Sort by max norm, strand, seq id, hom name and score */
   qsort(temp, tot, sizeof(struct hom), norm_cmp);

   /*--Initialize----------------------------------------------------*/
   /*-----------Rank counts results that are either potential TPs or
	TNs, results from the block are skipped ----------------------*/
   nhit = ntp = nblock = nfrag = t = itn = itp = rank = 0;
   roc = (double) 0.0;

   /*----------Process first sequence -------------------------------*/
   printf("\n\nTop %d sequences sorted by maximum normalized score:", MaxHit);
   printf("\n(the normalized score is the Score divided by the\n");
   printf("99.5 score in the block if available, or by the\n");
   printf("smallest Score in the search results file for the block)");
   if (SeqType == DNA)
      printf("\n\nBlock     Rank Frame Score      Location(bp) Sequence");
   else
      printf("\n\nBlock     Rank Frame Score      Location(aa) Sequence");
   if (t < tot && !temp[t].tp)
   {
      nhit++; tn[itn++] = rank; rank++;
      if ( (LisSeq+FragSeq) > 0)
         roc += (double) ntp / (double) (LisSeq + FragSeq);
      if (SeqType == DNA)
      {
	locfirst = compute_loc(temp[t].frame, temp[t].offset+1, temp[t].seqlen);
	loclast = compute_loc(temp[t].frame, temp[t].offset + temp[t].width,
			      temp[t].seqlen);
      }
      else
      {
	locfirst = temp[t].offset+1; loclast = temp[t].offset+temp[t].width;
      }
      printf("\n%d.---------------------------------------", nhit);
      printf("-------------------------------------------");
      printf("\n%8s  %4d  %2d    %4d  %7ld-%7ld %s %s",
	 temp[t].ac, temp[t].rank+1, temp[t].frame, temp[t].norm,
	 locfirst, loclast,
         temp[t].seq_id, temp[t].title);
   } /* end of not in the true positive list */
   else
   {	/*  In the TP list */
      if (temp[t].tp == 2) { nblock++; }
      else { ntp++; tp[itp++] = rank; rank++; }
   }
   strcpy(save_id, temp[t].seq_id);
   save_norm = temp[t].max_norm;
   if (temp[t].frame < 0) save_strand = -1; else save_strand = 1;
   save_t = t;
   t++;
   /*-------------------Print the rest of the hits -----------------*/
   /*  New hit if different seq_id or if same seq_id but different strand */
   while(nhit <= MaxHit && t <= tot)
   {
      if (temp[t].frame < 0) strand = -1; else strand = 1;
      if ( strcmp(temp[t].seq_id, save_id) != 0 ||
          (strcmp(temp[t].seq_id, save_id) == 0 &&
           strand != save_strand) )    /*  new group */
      {
	 if (blocks != NULL && !temp[save_t].tp)     /* process previous group */
	 {
	     if (blocks != NULL) check_dat(save_t, t, temp, blocks);
             else printf("\nBlocks not found in database\n");
	 }
	 strcpy(save_id,temp[t].seq_id);
         if (temp[t].frame < 0) save_strand = -1; else save_strand = 1;
         save_t = t;
	 if (t < tot && !temp[t].tp)
         {
            nhit++; tn[itn++] = rank; rank++; 
            if ( (LisSeq+FragSeq) > 0)
                roc += (double) ntp / (double) (LisSeq + FragSeq);
            if (nhit <= MaxHit)
            {
               printf("\n%d.---------------------------------------", nhit);
               printf("-------------------------------------------");
            }
	 }   /* end of if not a TP */
         else
         {
            if (temp[t].tp == 2) { nblock++; }
            else { ntp++; tp[itp++] = rank; rank++; }
         }
      }
      if (nhit <= MaxHit && t < tot && !temp[t].tp)
      {
         if (SeqType == DNA)
         {
	   locfirst = compute_loc(temp[t].frame,
		 temp[t].offset+1, temp[t].seqlen);
	   loclast = compute_loc(temp[t].frame, temp[t].offset + temp[t].width,
			      temp[t].seqlen);
         }
         else
         {
	   locfirst = temp[t].offset+1; loclast = temp[t].offset+temp[t].width;
         }
	 printf("\n%8s  %4d  %2d    %4d  %7ld-%7ld %s %s",
	    temp[t].ac, temp[t].rank+1, temp[t].frame, temp[t].norm,
	    locfirst, loclast,
            temp[t].seq_id, temp[t].title);
      }  /*  end of if not in the true positive list */
      t++;
   }

/*-------------------------------------------------------------------*/
   /*-- For statistics calcs, notice that the "scores" in the tn[] and
       tp[] arrays are really ranks & so are sorted in reverse order
       compared with matodat/blastdat/fastodat results; lower value
       is "higher" score --*/
   if (nhit > 0) roc /= (double) nhit;
   else roc = 1.0;

   tn999 = (int) ( (double) (SWISS-LisSeq-FragSeq-BlkSeq) * 0.001);
   if (tn999 > nhit) tn999 = nhit - 1;
   if (tn999 < 0) tn999 = 0;
   tp999 = 0;
   for (i=0; i < ntp; i++)
      if (tp[i] <= tn[tn999]) tp999++;

   pearson = LisSeq + FragSeq - ntp;
   i = ntp - 1;
   while (pearson >= 0 && pearson < nhit &&
          i >= 0 && i < ntp &&
          tn[pearson] < tp[i])
   {  pearson++; i--; }
   /*  ntp=#tps found not in block, nblock=#tps found in block,
       nhit=#found not in .lis file => tns for testing purposes */
   if (Fdat != NULL)
     fprintf(Fdat, "%s %d %d %d %d %d %d %d %d %d %d %d %f\n",
        AC, NScore, MaxHit, LisSeq, BlkSeq, FragSeq, LisSeq+FragSeq,
        nblock, ntp, nhit, tp999, pearson, roc);

}   /* end of show_hom */
/*======================================================================*/
/*   Sort by strand, seq id, rank */
int id_cmp(t1,t2)
struct hom *t1, *t2;
{
   int strand1, strand2;

   strand1 = strand2 = 1;
   if (t1->frame < 0) strand1 = -1;
   if (t2->frame < 0) strand2 = -1;

   if (strand1 != strand2) return(strand1 - strand2);
   else if (strcmp(t1->seq_id, t2->seq_id) > 0)  return(1);
   else if (strcmp(t1->seq_id, t2->seq_id) < 0)  return(-1);
   else                       return(t1->rank -  t2->rank);
}  /* end of id_cmp */
/*======================================================================*/
/*   Sort by strand, seq id, normalized score (d)  */
int idnorm_cmp(t1,t2)
struct hom *t1, *t2;
{
   int strand1, strand2;

   strand1 = strand2 = 1;
   if (t1->frame < 0) strand1 = -1;
   if (t2->frame < 0) strand2 = -1;

   if (strand1 != strand2) return(strand1 - strand2);
   else if (strcmp(t1->seq_id, t2->seq_id) > 0)  return(1);
   else if (strcmp(t1->seq_id, t2->seq_id) < 0)  return(-1);
   else                       return(t2->norm -  t1->norm);
}  /* end of idnorm_cmp */
/*======================================================================*/
/*   Sort by min_rank, strand, seq id, hom name */
int rank_cmp(t1,t2)
struct hom *t1, *t2;
{
   int strand1, strand2;

   strand1 = strand2 = 1;
   if (t1->frame < 0) strand1 = -1;
   if (t2->frame < 0) strand2 = -1;

   if (t1->min_rank != t2->min_rank)  return(t1->min_rank - t2->min_rank);
   else if (strand1 != strand2)  return(strand1 - strand2);
   else if (strcmp(t1->seq_id, t2->seq_id) > 0)          return(1);
   else if (strcmp(t1->seq_id, t2->seq_id) < 0)          return(-1);
   else                               return(strcmp(t1->ac, t2->ac));
}  /* end of rank_cmp */
/*======================================================================*/
/*   Sort by max normalized score (d), strand, seq id, hom name */
int norm_cmp(t1,t2)
struct hom *t1, *t2;
{
   int strand1, strand2;

   strand1 = strand2 = 1;
   if (t1->frame < 0) strand1 = -1;
   if (t2->frame < 0) strand2 = -1;

   if (t1->max_norm != t2->max_norm)  return(t2->max_norm - t1->max_norm);
   else if (strand1 != strand2)  return(strand1 - strand2);
   else if (strcmp(t1->seq_id, t2->seq_id) > 0)          return(1);
   else if (strcmp(t1->seq_id, t2->seq_id) < 0)          return(-1);
   else                               return(strcmp(t1->ac, t2->ac));
}  /* end of norm_cmp */
/*======================================================================*/
void print_blurb()
{
   FILE *fstp;
   char line[MAXLINE];

   if ((fstp=fopen("blksrch.stp", "r")) == NULL)
   {
      printf("\n========================================");
      printf("=======================================");
      printf("\nSearch results from the BLOCKS e-mail searcher.");
      printf("\nPlease report problems to jorja@sparky.fhcrc.org");
      printf(", include your query\nand this output.");
      printf(" To obtain help, send the word HELP on a single");
      printf("\nline to blocks@howard.fhcrc.org");
      printf("\n========================================");
      printf("=======================================");
      printf("\nCopyright (c) 1992 by the Fred Hutchinson Cancer Research Center");
      printf("\nIf you use BLOCKS in your research, please cite:");
      printf("\nSteven Henikoff and Jorja G. Henikoff,");
      printf(" Automated assembly of protein");
      printf("\nblocks for database searching, NAR 19:23 (1991), 6565-6572.");
/*   printf("\n2. James C. Wallace and Steven Henikoff,");
      printf(" PATMAT: a searching and extraction");
      printf("\n   program for sequence, pattern and block queries");
      printf(" and databases,");
      printf("\n   CABIOS 8:3 (1992), 249-254.");
*/
      printf("\n========================================");
      printf("=======================================");
      printf("\nEach numbered result consists of one");
      printf(" or more blocks from a PROSITE group");
      printf("\nfound in the query sequence. One set");
      printf(" of the highest-scoring blocks that");
      printf("\nare in the correct order and separated");
      printf(" by distances comparable to the BLOCKS");
      printf("\ndatabase is selected for analysis.");
      printf(" If this set includes multiple blocks");
      printf("\nthe probability that the lower scoring");
      printf(" blocks support the highest scoring");
      printf("\nblock is reported. Maps of the database");
      printf(" blocks and query sequence are shown:");
      printf("\n  < indicates the sequence has been");
      printf(" truncated to fit the page");
      printf("\n  : indicates the minimum distance");
      printf(" between blocks in the database");
      printf("\n  . indicates the maximum distance");
      printf(" between blocks in the database");
      printf("\nThe maps are aligned on the highest");
      printf(" scoring block. The alignment of the");
      printf("\nquery sequence with the sequence");
      printf(" closest to it in the BLOCKS database");
      printf("\nis shown. Upper case in the query");
      printf(" sequence indicates at least one");
      printf("\noccurrence of the residue in that");
      printf(" column of the block.");
      printf("\n========================================");
      printf("=======================================");
      printf("\n");
   }
   else
   {
      while (!feof(fstp) && fgets(line, MAXLINE, fstp) != NULL)
         printf("%s", line);
      fclose(fstp);
   }

}  /* end of print_blurb */
/*======================================================================*/
/*   Look up this group in blocks.dat                                   */
/*======================================================================*/
void check_dat(min_t, max_t, results, fblock)
int min_t, max_t;
struct hom *results;
struct best_block *fblock;
{
   int nblock;

   /*---- Add the query as the last sequence in each block in the list ---*/
   nblock = add_query(min_t, max_t, results, fblock);

   /*--- Map the blocks if more than one in the hit ----*/
   if (nblock > 1) map_blocks(min_t, max_t, results, fblock);
   else            printf("\n");

   /*---Determine the closest block seq to the query seq & display them --*/
   align_blocks(fblock);

}  /*  End of check_dat */
/*=====================================================================
      Read the blocks database
      Truncates blocks wider than MAX_WIDTH
=======================================================================*/
struct best_block *read_block(fblk, ac)
FILE *fblk;
char *ac;
{
   char line[MAXLINE], *ptr, *ptr1;
   struct best_block *block;
   int done, aclen;

   done = NO;
   block = NULL;
   aclen = strlen(ac);
   while (!done && !feof(fblk) && fgets(line, MAXLINE, fblk) != NULL)
   {
      if (strncmp(line, "AC   ", 5) == 0 &&
          strncmp(line+5, ac, aclen) > 0)
	  done = YES;
      else if (strncmp(line, "AC   ", 5) == 0 &&
               strncmp(line+5, ac, aclen) == 0)
      {
	 block = (struct best_block *) malloc(sizeof(struct best_block));
	 if (block == NULL)
	 {  printf("\nOUT OF MEMORY\n\n");  exit(-1);  }
         block->query = block->anchor = NO;
	 block->next_block = NULL;
	 block->dat_pos = ftell(fblk);
	 strncpy(block->ac, line+5, aclen + 1); block->ac[aclen + 1] = '\0';
	 if (block->ac[aclen] == ';') block->ac[aclen] = '\0';
	 block->minprev = block->maxprev = -1;
	 block->query = NO; block->rank = 9999;
	 ptr = strtok(line+12, "(");
	 if (ptr != NULL)
	 {
	    ptr = strtok(NULL, ",");
	    if (ptr != NULL)
	    {
	       block->minprev = atoi(ptr);
	       ptr = strtok(NULL, ")");
	       if (ptr != NULL) block->maxprev = atoi(ptr);
	    }
	 }
      }
      else if (block != NULL && strncmp(line, "BL   ", 5) == 0)
      {
	 block->s995 = 0;
	 ptr=strstr(line,"99.5\%=");
	 if (ptr != NULL)
	 {
	    ptr1 = strtok(ptr, "="); ptr1=strtok(NULL, "\n\r");
	    block->s995 = atoi(ptr1);
	 }
	 block->strength = 0;
	 ptr=strstr(line,"strength=");
	 if (ptr != NULL)
	 {
	    ptr1 = strtok(ptr, "="); ptr1=strtok(NULL, "\n\r");
	    block->strength = atoi(ptr1);
	 }
	 fill_block(fblk, block);
	 return(block);
      }
   }
   return(NULL);		/* no block found! */
}  /* end of read_block */
/*====================================================================
     Fill up the block structure
=======================================================================*/
void fill_block(fblk, block)
FILE *fblk;
struct best_block *block;
{
   int done, i, n, cluster, ncluster;
   char line[MAXLINE], ctemp[MAXLINE], *ptr, *ptr1;

   block->nseq = block->width = 0;
   cluster = ncluster = 0;
   done=NO;
   while (!done && !feof(fblk) && fgets(line,MAXLINE,fblk) != NULL)
   {
      if (strlen(line) == 1)		/* blank line => new cluster */
      {
	 /*  Set #seqs in cluster to seqs in previous cluster */
	 if (ncluster > 0)
	   for (n=0; n<block->nseq; n++)
	      if (block->cluster[n] == cluster) block->ncluster[n] = ncluster;
	 cluster++; ncluster = 0;
      }
      else if ((int) strlen(line) > (int) 1)
      {
	 if (strncmp(line, "//", 2) == 0) done=YES;
	 else if ((int) strlen(line) > 20)
	 {
	    ptr=strtok(line, "(");
	    /*  need to strip leading spaces off ptr here  */
	    strcpy(ctemp, ptr);
	    i=0;
	    while (ctemp[i] == ' ') i++;
	    strcpy(block->name[block->nseq], ctemp+i);
	    ptr=strtok(NULL, ")");
	    block->offset[block->nseq] = atoi(ptr);

	    ptr1=strtok(NULL, "\n\r");
	    i=0;
	    while (ptr1[i] == ' ') i++;
	    ptr = strtok(ptr1+i, " \t\r\n");

            /* Truncate wide blocks   */
            if (strlen(ptr) > MAX_WIDTH) ptr[MAX_WIDTH - 1] = '\0';
	    strcpy(block->aa[block->nseq], ptr);
	    block->cluster[block->nseq] = cluster;
	    ncluster++;		/* # seqs in current cluster */
	    block->width = (int) strlen(block->aa[block->nseq]);
	    block->nseq++;
	 }
      }
   }
   /*  Compute weights for the last cluster */
   if (ncluster > 0)
	   for (n=0; n<block->nseq; n++)
	      if (block->cluster[n] == cluster) block->ncluster[n] = ncluster;
}  /* end of fill_block */
/*=======================================================================
     Add the query sequence as the last sequence in each block for
     which it has a hit.
==========================================================================*/
int add_query(min_t, max_t, results, fblock)
int min_t, max_t;
struct hom *results;
struct best_block *fblock;
{
   struct best_block *block, *maxblock;
   int i, t, j, imin, nblock, mblock, lastrank, one;
   char ac[10], pline[MAXLINE], tline[40];
   struct temp *temp;
   double prob;

   /*------- Need to know how many possible blocks there are ---*/
   /*--------Clear out any previous sequence information -------*/
   nblock=0;
   block = fblock;
   while (block != NULL)
   {
      nblock++;
      block->name[block->nseq][0] = '\0';
      block->offset[block->nseq] = 0;
      block->cluster[block->nseq] = -1;
      block->ncluster[block->nseq] = -1;
      block->aa[block->nseq][0] = '\0';
      if (block->query) block->nseq--;
      block->query = NO;
      block = block->next_block;
   }
   if (nblock == 1) one = YES;
   else             one = NO;

   /*------ Sort this group of hits by rank now --------------*/
   temp = (struct temp *) malloc((max_t-min_t)*sizeof(struct temp));
   if (temp == NULL)
   {
      printf("\nOUT OF MEMORY\n\n"); exit(-1);
   }
   for (i=0; i<max_t-min_t; i++)
   {
      temp[i].value = results[i+min_t].rank;
      temp[i].index = i+min_t;	temp[i].flag = 0;
   }
   qsort(temp, max_t-min_t, sizeof(struct temp), tempcmp);

   /*----------- Piece together a set of compatible hits, highest
       rank first and compute probability of multiple blocks ----------*/
   prob = 1.0; maxblock = NULL; pline[0] = '\0';
   mblock = 0;  	/* number of blocks in the map */
   lastrank = 0;	/* rank of previous block in the map */
   for (i=0; i<max_t-min_t; i++)
   {
      t = temp[i].index;
      block = fblock;
      while (block != NULL && t < max_t)
      {
         DBLen = results[t].seqlen;
	 strcpy(ac, results[t].ac);
	 if (strcmp(ac, block->ac) == 0 && !block->query &&
             (i==0 ||
	     distance_okay(results[t].offset, fblock, block)) )
	 {
	    /*--- Adds the query to the block  ------*/
	    strncpy(block->name[block->nseq], results[t].seq_id, SNAMELEN);
            block->name[block->nseq][SNAMELEN] = '\0';
	    block->offset[block->nseq] = results[t].offset;
	    block->cluster[block->nseq] = -1;
	    block->ncluster[block->nseq] = 1;
            /*   Search alignment may be shorter than block ... 
                 pad with spaces */
            if ((int) strlen(results[t].aa) < block->width )
            {
               for (j = strlen(results[t].aa); j < block->width; j++)
                      results[t].aa[j] = ' ';
               results[t].aa[block->width] = '\0';
            }
	    strcpy(block->aa[block->nseq], results[t].aa);
               
	    block->query = results[t].score;
	    block->rank = results[t].rank;
	    block->nseq++;
            results[t].map_flag = YES;
	    /*--- Compute probability -----*/
            if (i==0)
            {
                maxblock = block;
                maxblock->anchor = YES;
                lastrank = temp[i].value;
            }
            else if (i>0 && NRank > 0.0 && DBLen > 0.0)
            {
               mblock++;     /* number of supporting blocks so far */
               block->anchor = NO;

/* blksort calcs: compute prob a block of this family could rank this high
               n = (double) DBType*DBLen*NBlock-lastrank;
	       r = (double) temp[i].value-lastrank;
	       p = (double) nblock-mblock;
	       prob *= hypergeo(n, r, p, 1.0);
*/
               /*   Probability supporting block ranks this high */
               prob *= (double) (temp[i].value + 1) / NRank;
               /* Compute prob this block is this far from the anchor block */
               /* DBLen is always in protein units, distance is always
                  in protein units no need to multiply by DBType to convert to
                  correct units. Min dist = -1; Max dist = distance() */
               prob *= (double) (distance(maxblock,block)+1) / DBLen;
               strcat(pline, block->ac); strcat(pline, " ");
               lastrank = temp[i].value;
            }
	 }
/* blksort calcs         else if (strcmp(ac, block->ac) == 0)
            nblock--;                           block doesn't fit 
*/
	 block = block->next_block;
      }
   }
   if (prob > 1.0) prob = 1.0;
   if (prob < 0.0) prob = 0.0;
   if (maxblock != NULL && strlen(pline))
   {
      printf("\n\nP<%6.2g for ", prob);
      i=0;
      while (i < (int) strlen(pline))
      {
         imin=36;
         if (i+imin > (int) strlen(pline)) imin = (int) strlen(pline) - i;
         strncpy(tline, pline+i, imin); tline[imin] = '\0';
         if (i==0) printf("%s", tline);
         else printf("\n                %s", tline);
         i += 36;
      }
      printf("in support of %s", maxblock->ac);
   }
   else if (maxblock != NULL)
   {   printf("\n\nNo P-value computed for single block hits."); }
   else
   {   printf("\n\nERROR: Anchor block not found in blocks database."); }

   /*----- multimat.dat statistics ------------------------------*/
   for (i=0; i<max_t-min_t; i++)
   {
      t = temp[i].index;
      if (results[t].map_flag && Fout != NULL)
         fprintf(Fout, "%s %s %d %d %d %d %d %d %6.2g %s\n",
           results[t].ac, results[t].seq_id, results[t].rank + 1,
           results[t].frame, results[t].score, results[t].offset + 1,
           nblock, mblock+1, prob, results[t].aa);
   }

   free(temp);
   return(mblock+1);
}  /*  end of add_query */
/*===================================================================
    Compute probability from hypergeometric distribution:
    Assume there are n objects, p of one kind and n-p of another.
    Then when r objects are selected at random from among the n
    objects, the prob. that k of the r are of type p is:
  
      ( p )  ( n-p )        p!               (n-p)!
      ( k )  ( r-k )       ----- * -----------------------
                          k!(p-k)!  (r-k)!( (n-p)-(r-k) )!
    -----------------  =  ---------------------------------
         ( n )                       n!
         ( r )                   ----------
                                  r!(n-r)!

       p!         r!          (n-p)!         (n-r)!
  =  -------- * --------  * ----------  * -----------------
     k!(p-k)!    (r-k)!         n!         ( (n-r)-(p-k) )!

  Which can be simplified to a product of these expressions:
   1.  k terms:
       p*(p-1)*(p-2)*...(p-k+1) / k!
   2.  k terms:
       r*(r-1)*(r-2)*...(r-k+1)
   3.  p terms:
       1/(n*(n-1)*(n-2)*...(n-p+1))
   4.  p-k terms:
       (n-r)*(n-r-1)*....(n-r-(p-k)+1)
      
=====================================================================*/
double hypergeo(n, r, p, k)
double n, r, p, k;
{
   double prob;
   int i;

   prob = 1.0;
   for (i=0; i<k; i++)  prob *= (double) (p-i) * (r-i) / (i+1);
/*   for (i=0; i< (p-k); i++) prob *= (double) (n-r-i);
   for (i=0; i<p; i++) prob /= (double) (n-i);
*/
   for (i=0; i<(p-k); i++) prob *= (double) (n-r-i)/(n-i);
   for (i=(p-k); i<p; i++) prob /= (double) (n-i);
 
   return(prob);
}  /* end of hypergeo */
/*====================================================================
     Check if this block will fit between existing blocks
       Minimum distance between blocks is assumed to be 0 and
       maximum distance is minprev+maxprev = 2xthe average distance
       between the blocks in sequences in the blocks database.
=====================================================================*/
int distance_okay(offset, fblock, nblock)
int offset;
struct best_block *fblock, *nblock;
{
   struct best_block *block, *lblock, *rblock;
   int mindist, maxdist;

   lblock = rblock = NULL;
   block = fblock;
   while (block != NULL && rblock == NULL)
   {
      if (block->query)		/* query added to this block already */
      {
	 if (strcmp(block->ac, nblock->ac) < 0)
	    lblock = block;	/* closest block on the left */
	 else if (strcmp(block->ac, nblock->ac) > 0 && rblock == NULL)
	    rblock = block;	/* closest block on the right */
      }
      block = block->next_block;
   }
   /*---------  Check the distance to the left --------------------*/
   if (lblock != NULL)
   {
      mindist = maxdist = lblock->offset[lblock->nseq-1] + lblock->width;
      block = lblock->next_block;
      while (block != NULL && block != nblock)
      {
	 if (strcmp(block->ac, nblock->ac) < 0)
         {
            mindist += block->width;
	    maxdist += prev_dist(block) + block->width;
         }
	 block = block->next_block;
      }
      if (mindist-2 > offset) return(NO); /* allow overlap of 2 */
      if (maxdist + prev_dist(nblock) < offset) return(NO);
   }
   /*----------Check the distance to the qright ---------------------*/
   if (rblock != NULL)
   {
      mindist = maxdist = offset + nblock->width;
      block = nblock->next_block;
      while (block != NULL && block != rblock)
      {
	 if (strcmp(block->ac, rblock->ac) < 0)
         {
            mindist += block->width;
	    maxdist += prev_dist(block) + block->width;
         }
	 block = block->next_block;
      }
      if (mindist-2 > rblock->offset[rblock->nseq-1]) return(NO);
      if (maxdist + prev_dist(rblock) < rblock->offset[rblock->nseq-1])
	 return(NO);
   }
   return(YES);
}  /*  end of distance_okay */
/*===================================================================
    Compute maximum distance preceding a block
=====================================================================*/
int prev_dist(block)
struct best_block *block;
{
   int dist;

/*
   dist = block->minprev + block->maxprev;
*/
   dist = 2 * block->maxprev;
   if (dist < 1) dist = 1;
   return(dist);
}  /* end of prev_dist */
/*===================================================================
     Compute maximum allowable distance between two blocks
====================================================================*/
int distance(fromblock, toblock)
struct best_block *fromblock, *toblock;
{
    int dist, maxdist;
    struct best_block *block, *lblock, *rblock;
  
    if (fromblock == NULL || toblock == NULL) return(-1);

    if (strcmp(fromblock->ac, toblock->ac) < 0)
    { lblock=fromblock; rblock=toblock;
      maxdist = DBLen-(fromblock->offset[fromblock->nseq-1]+fromblock->width);
    }
    else
    { lblock=toblock; rblock=fromblock;
      maxdist = fromblock->offset[fromblock->nseq-1];
    }
/*     For BLIMPS matrix search, DBLen is always in protein units
    if (DBType*maxdist > DBLen) maxdist = (int) DBLen/DBType;
*/
    if (maxdist > DBLen) maxdist = (int) DBLen;

    dist = 0;
    block=lblock->next_block;
    while (block != NULL && block != rblock)
    {
       dist += prev_dist(block) + block->width;
       block = block->next_block;
    }
    dist += prev_dist(rblock);
    /*  Don't want possible distance to be larger than possible
        sequence positions! */
    if (dist > maxdist) dist = maxdist;
    return(dist);
}  /* end of distance */
/*====================================================================
    Map a path of blocks
    The mapping is done in three phases:
    1. The blocks in the database are mapped from the 1st position of
       the first block to the last position of the last block. This
       range divided by MAXMAP (number of columns for the map to occupy
       when printed) determines the scale of the map.
    2. The blocks in the query sequence that are consistent with the
       database are mapped to the same scale as the database blocks.
       The two maps are aligned on the highest scoring block (maxblock).
       The query map may occupy more than MAXMAP columns and may have
       to be truncated or compressed. A maximum of MAXLINE characters
       are mapped before compression.
    3. The blocks in the query sequence that are not consistent with the
       database are mapped to the same scale. Since these blocks may
       overlap one another, they may take more than one line. Each line
       is aligned with the left end of the map in 2.
=======================================================================*/
void map_blocks(min_t, max_t, results, fblock)
int min_t, max_t;
struct hom *results;
struct best_block *fblock;
{
   int totblks, maxdist, i, imin, imax, pos, maxspot, firstpos, lastpos;
   int spot, qspot, firstspot, lastspot, t, done, qleft, qright;
   double scale;
   struct best_block *block, *maxblock;
   char ctemp[20], dbline[MAXLINE], qline[MAXLINE], pline[MAXLINE];

   totblks = maxdist = 0;
   maxblock = NULL;
   /*--- Determine the scale for the database blocks (# aas per space) ---*/
   block = fblock;
   while (block != NULL)
   {
      totblks++;
      if (block != fblock)
         maxdist += block->maxprev;  /* distance from last block */
      maxdist += block->width;    /* width of this block */
      if (block->query && block->anchor) maxblock = block;
      block = block->next_block;
   }
   scale = (double) maxdist/MAXMAP;       /* max of 60 spaces per line */

   for (spot=0; spot<MAXLINE; spot++)
   {  dbline[spot] = qline[spot] = ' '; }

   /*---  1. Now map the database blocks ------------------------------*/
   spot=maxspot=0;
   block = fblock;
   while (block != NULL)
   {
      if (block != fblock)
      {
         imin = (int) ((double) 0.5+block->minprev/scale);
         imax = (int) ((double) 0.5+(block->maxprev-block->minprev)/scale);
         for (i=spot; i < spot+imin; i++) dbline[i] = ':';
         spot += imin;
         for (i=spot; i < spot+imax; i++) dbline[i] = '.';
         spot += imax;
      }
      imin = (int) ((double) 0.5+block->width/scale);
      if (imin < 1) imin = 1;
      for (i=spot; i < spot+imin; i++) strncpy(dbline+i, block->ac+ACLen, 1);
      spot += imin;
      if (block==maxblock) maxspot=spot;
      block = block->next_block;
   }
   dbline[spot] = '\0';
   
   printf("\n                         |-----%5d residues----|",
          (int) ((double) 0.5+(scale*25.0)) );
   printf("\n%-20s ", AC);		/* need SNAMELEN here to line up */
   printf("%s", dbline);

   /*---- 2. Now map the query sequence if there is one ----------------*/
   /*---Line up the highest scoring query block with the database 
        pos = sequence position,
        spot=scaled pos in qline , qspot=scaled pos of anchor block in query,
        maxspot=scaled pos of anchor block in block map */
   /* firstpos & lastpos are first & last sequence positions mapped into qline
      firstspot is position of firstpos in qline, etc. */
   spot=qspot=pos=0;
   firstpos = lastpos = firstspot = lastspot = -1;
   block = fblock;
   while (block != NULL)
   {
      if (block->query)
      {
         if (firstpos < 0) firstpos = block->offset[block->nseq-1];
         lastpos = block->offset[block->nseq-1] + block->width;
	 imin = (int) ((double) 0.5+(block->offset[block->nseq-1]-pos)/scale);
         if (imin >= 0)   /* skip over blocks in wrong order */
         {
	    if (spot+imin > MAXLINE && spot+1 < MAXLINE)
               qline[spot++] = '<';
            else
            {
	       for (i=spot; i<spot+imin; i++) qline[i] = ':';
	       spot += imin;
            }
	    imin = (int) ((double) 0.5+block->width/scale);
	    if (imin < 1) imin = 1;
            if (spot+imin < MAXLINE)
            {
	       for (i=spot; i<spot+imin; i++) 
                   strncpy(qline+i, block->ac+ACLen, 1);
	       spot += imin;
            }
	    pos = block->offset[block->nseq-1] + block->width;
	    if (block == maxblock) qspot = spot;
         }
      }
      block = block->next_block;
   }
   qline[spot] = '\0';

   /*  Line up qline/qspot with dbline/maxspot */
   strncpy(ctemp, results[min_t].seq_id, SNAMELEN); ctemp[SNAMELEN] = '\0';
   printf("\n%-20s ", ctemp);
   if (qspot <= maxspot)
   {
       qleft = 0; firstspot = maxspot - qspot;
       for (i=qspot; i<maxspot; i++) printf(" ");
   }
   else
   {
      qleft = qspot - maxspot; firstspot = 0;
      qline[0] = '<';
   }
   if (((int) strlen(qline) - qspot) <= ((int) strlen(dbline) - maxspot))
      qright = (int) strlen(qline);
   else
      qright = qspot + (int) strlen(dbline) - maxspot;
   strncpy(pline, qline+qleft, qright-qleft+1); pline[qright-qleft+1] = '\0';
   if (qleft > 0 && (pline[0] == '.' || pline[0] == ':')) pline[0] = '<';
   if (qright != (int) strlen(qline)) pline[(int) strlen(pline)-1] = '>';
   printf("%s\n", pline);

   /*--------3. Now map the inconsistent hom hits wrt consistent ones--*/
   /* Take the hits as they come, as many as will fit on a line        */
   done = NO; 
   while (!done)
   {
      done=YES; spot=0;
      for (i=0; i<MAXLINE; i++) qline[i] = ' ';
      lastspot = 0;
      for (t=min_t; t<max_t; t++)
      {
         if (!results[t].map_flag)
         {
            imin = (int) ((double) 0.5+(results[t].offset)/scale);
            if (results[t].offset < firstpos) spot = firstspot - imin;
            else spot = firstspot + imin;
            if (spot >= 0 && spot < MAXLINE && qline[spot] == ' ')
            {
               imin = (int) ((double) 0.5+results[t].width/scale);
               if (imin < 1) imin = 1;
               i = spot;
               while (i < spot+imin && i < MAXLINE)
               {
                  strncpy(qline+i, results[t].ac+ACLen, 1);
                  if (i > lastspot) lastspot = i;
	          i++;
               }
               done = NO;
               results[t].map_flag = YES;
            }
         }
      }  /* end of a line of output */
      qline[lastspot+1] = '\0';
      if (!done && strlen(qline))
      {
         printf("%-20s %s\n", ctemp, qline); /* ctemp is hit name */
      }
   } /* end of inconsistent hits */

   printf("\n");

}  /* end of map_blocks */
/*========================================================================
	Align the query sequence with the sequence closest to it in
	each block for which it has a hit.
	datline from blocks database
	homline from blimps search results
=========================================================================*/
void align_blocks(fblock)
struct best_block *fblock;
{
   struct best_block *block;
   int s, i, itemp, spot, bspot, datmin, datmax, homdist;
   char datline[MAXLINE], homline[MAXLINE], blkline[MAXLINE];
   char barline[MAXLINE], saveac[MAXAC+1];
   char ctemp[MAXLINE];

   for (i=0; i<MAXLINE; i++)
   { datline[i] = homline[i] = blkline[i] = barline[i] = ' ';}
   bspot=spot=datmin=datmax=homdist=0;
   saveac[0] = '\0';
   for (i=0; i<ACLen+1; i++) strcat(saveac, " ");
   block=fblock;
   while (block != NULL)
   {
      if (block->query)
      {
	 datmin += block->minprev;
	 datmax += block->maxprev;
	 homdist = block->offset[block->nseq-1] - homdist;
         consensus(block);
/*>>> if aborts on call to closest_seq(), reduce MAXSEQS in motifj.h
	or set s = 0
<<<*/
	 s = closest_seq(block);
	 bspot = block->width;
	 /*--- 28 spaces for name & offset on datline & homline;
	       27 spaces for distance info on blkline ----*/
	 if (bspot < 27) bspot = 27;   /* min blkline width */
	 if ((spot + bspot + 28) > 80) 
         {
            itemp = spot;
            if (bspot > itemp) itemp = bspot;
	    homline[itemp]= datline[itemp]= blkline[itemp]= barline[itemp]= '\0';
	    printf("\n%s", blkline);
            printf("\n%s", datline); printf("\n%s", barline);
	    printf("\n%s\n", homline);
            for (i=0; i<MAXLINE; i++)
	    { datline[i] = homline[i] = blkline[i] = barline[i] = ' '; }
            spot=0;
         }
	 sprintf(ctemp, "%-20s", block->ac);		/* SNAMELEN */
	 strncpy(blkline+spot, ctemp, strlen(ctemp));
         sprintf(ctemp, "%-20s", block->name[s]);	/* SNAMELEN */
	 strncpy(datline+spot, ctemp, strlen(ctemp));
	 sprintf(ctemp, "%-20s", block->name[block->nseq-1]);
	 strncpy(homline+spot, ctemp, strlen(ctemp));
	 spot+=SNAMELEN + 1;   /* name width + 1 */
/*
         strncpy(blkline+spot, block->ac, strlen(block->ac));
         strncpy(datline+spot, block->name[s], strlen(block->name[s]));
         strncpy(homline+spot, block->name[block->nseq-1],
              strlen(block->name[block->nseq-1]));
*/

	 sprintf(ctemp, "%c<->%c", saveac[ACLen], block->ac[ACLen]);
	 strncpy(blkline+spot, ctemp, strlen(ctemp));
	 sprintf(ctemp, "%6d", block->offset[s]);
	 strncpy(datline+spot, ctemp, strlen(ctemp));
	 sprintf(ctemp, "%6d", block->offset[block->nseq-1] + 1);
	 strncpy(homline+spot, ctemp, strlen(ctemp));
	 spot+=9;   /* offset width + 2 */
/*
	 blkline[spot] = saveac[ACLen];
	 strncpy(blkline+spot+1, "<->", 3);
	 blkline[spot+4] = block->ac[ACLen];
	 kr_itoa(block->offset[s], ctemp, 10);
	 strncpy(datline+spot, ctemp, strlen(ctemp));
	 kr_itoa(block->offset[block->nseq-1]+1, ctemp,10);
	 strncpy(homline+spot, ctemp, strlen(ctemp));
*/

	 bspot = spot;
         sprintf(ctemp, "(%d,%d):%d", datmin, datmax, homdist);
/*
	 strncpy(blkline+bspot, "(", 1); bspot++;
	 kr_itoa(datmin, ctemp, 10);
	 strncpy(blkline+bspot, ctemp, strlen(ctemp));
	 bspot += strlen(ctemp);
	 strncpy(blkline+bspot, ",", 1); bspot++;
	 kr_itoa(datmax, ctemp, 10);
	 strncpy(blkline+bspot, ctemp, strlen(ctemp));
	 bspot += strlen(ctemp);
	 strncpy(blkline+bspot, ")", 1); bspot++;
	 strncpy(blkline+bspot, ":", 1); bspot++;
	 kr_itoa(homdist, ctemp, 10);
*/
	 strncpy(blkline+bspot, ctemp, strlen(ctemp));
	 bspot += strlen(ctemp);

	 for (i=0; i<block->width; i++)
         {
	    strncpy(datline+spot, block->aa[s]+i, 1);
	    strncpy(homline+spot, block->aa[block->nseq-1]+i, 1);
            if (strncmp(datline+spot, homline+spot, 1) == 0)
                barline[spot] = '|';
	    spot++;
         }
	 spot+=3;	/* space between blocks on same line */
	 datmin = datmax = 0;
	 homdist = block->offset[block->nseq-1] + block->width;
	 strcpy(saveac, block->ac);
      }
      else
      {		/*  block not in query sequence */
	 datmin += block->minprev + block->width;
	 datmax += block->maxprev + block->width;
      }
      block = block->next_block;
   }
   itemp = spot;
   if (bspot > itemp) itemp = bspot;
   homline[itemp] = datline[itemp] = blkline[itemp] = barline[itemp] = '\0';
   printf("\n%s", blkline);
   printf("\n%s", datline); printf("\n%s", barline);
   printf("\n%s\n", homline);
}  /* end of align_blocks */
/*======================================================================*/
/*--- Just find the sequence closest to the last sequence in the block
      based on the number of identities */
/*======================================================================*/
int closest_seq(block)
struct best_block *block;
{
   int npair, s1, s2, l1, l2, i, i1, i2;
   int maxscore, maxs1;
   struct pair pairs[MAXSEQS*(MAXSEQS-1)/2]; 

   npair = block->nseq;

/*    Compute scores for all possible pairs of sequences            */
   for (s1=0; s1<block->nseq-1; s1++)   		/* col = 0, n-2     */
   {
      l1 = 0;
      s2 = block->nseq-1;   /* last seq in block == the query sequence */
      l2 = 0;
      pairs[s1].score = 0;
      pairs[s1].cluster = -1;
      for (i=0; i<=block->width; i++)
      {
	 i1 = l1+i;  i2 = l2+i;
	 /* s1 is the block seq & is all caps, but s2 is the query
	    seq & it may have lower case chars in it */
	 if (i1 >= 0 && i1 < block->width &&
	     i2 >= 0 && i2 < block->width &&
		strncasecmp(block->aa[s1]+i1, block->aa[s2]+i2, 1) == 0)
		   pairs[s1].score += 1;
      }
   }  /* end of s1 */
   maxs1 = -1; maxscore = -1;
   for (s1=0; s1 < block->nseq-1; s1++)
      if (pairs[s1].score > maxscore)
      {
	 maxscore = pairs[s1].score;
	 maxs1 = s1;
      }

      return(maxs1);
}  /*  end of cluster_seqs */
/*====================================================================*/
int tempcmp(t1, t2)
struct temp *t1, *t2;
{
   return(t1->value - t2->value);
}  /* end of tempcmp */
/*====================================================================*/
struct best_block *get_blocks(fblk, ac)
FILE *fblk;
char *ac;
{
   struct best_block *block, *fblock, *lblock;
   int nblk;

   nblk = 0;
   block = fblock = lblock = NULL;

   /*----- Make a list of blocks in the database -----*/
   while ( (block = read_block(fblk, ac)) != NULL)
   {
      if (nblk > 0) lblock->next_block = block;
      else fblock = block;
      nblk++;
      lblock=block;
   }
   if (fblock == NULL)
      printf("\nNo blocks for family %s found in database.", ac);
   return(fblock);
}  /* end of get_blocks */
/*======================================================================
    If the search was of a matrix database, the sequence segment will be
    all lower case. Change any residue that matches any sequence in the
    block to upper case. The block segments are all upper case.
========================================================================*/
void consensus(block)
struct best_block *block;
{
   int s, s1, i;
   char seqaa[2];

   s = block->nseq - 1;		/* This is the query sequence */
   for (s1=0; s1 < s; s1++)
      for (i=0; i < block->width; i++)
      {
         if (strncasecmp(block->aa[s]+i, block->aa[s1]+i, 1) == 0)
         {
            strncpy(seqaa, block->aa[s]+i, 1);
            if (seqaa[0] >= 97 && seqaa[0] <= 122)
            {
               seqaa[0] -= 32;
               strncpy(block->aa[s]+i, seqaa, 1);
            }
         }
      }

}  /* end of consensus */
/*=========================================================================
       Blimps reports seqlen in AAs for Block vs DNA search

  frame  1:1 => bp 1,     2:1 => 2,         3:1 => 3
           bp = frame + 3 * offset
       frame -1:1 => bp QLen, -2:1 => QLen - 1, -3:l => QLen - 2 
           bp = frame + (QLen + 1) - 3 * offset
=========================================================================*/
int compute_loc(frame, offset, seqlen)
int frame;
int offset;
double seqlen;
{
   int loc;

   loc = (offset - 1) * SeqType;
   if (frame < 0) loc = (seqlen * SeqType + 1) - loc;
   loc += frame;
   if (frame == 0) loc++;	/* Blimps frame = 0 for protein */
   return(loc);
}   /* end of compute_loc */
/*========================================================================
       If the blocks have 99.5 score in them, renormalize with it
=========================================================================*/
void re_normalize(temp, totscores, blocks)
struct hom *temp;
int totscores;
struct best_block *blocks;
{
   int i;
   struct best_block *block;

   for (i=0; i < totscores; i++)
   {
      block = blocks;
      while (block != NULL)
      {
         if (strcmp(temp[i].ac, block->ac) == 0 &&
             block->s995 > 0)
         {
      temp[i].norm = round((float) 1000. * temp[i].score / block->s995);
      temp[i].max_norm = temp[i].norm;
         }
         block = block->next_block;
      }
   }
}   /* end of re_normalize */