comparison clustalomega/clustal-omega-1.0.2/src/squid/sqio.c @ 1:bc707542e5de

Uploaded

1:bc707542e5de

/*****************************************************************
 * SQUID - a library of functions for biological sequence analysis
 * Copyright (C) 1992-2002 Washington University School of Medicine
 * 
 *     This source code is freely distributed under the terms of the
 *     GNU General Public License. See the files COPYRIGHT and LICENSE
 *     for details.
 *****************************************************************/

/* File: sqio.c
 * From: ureadseq.c in Don Gilbert's sequence i/o package
 *
 * Reads and writes nucleic/protein sequence in various
 * formats. Data files may have multiple sequences.
 *
 * Heavily modified from READSEQ package
 * Copyright (C) 1990 by D.G. Gilbert
 * Biology Dept., Indiana University, Bloomington, IN 47405
 * email: [e-mail]
 * Thanks Don!
 *
 * SRE: Modifications as noted. Fri Jul  3 09:44:54 1992
 *      Packaged for squid, Thu Oct  1 10:07:11 1992
 *      ANSI conversion in full swing, Mon Jul 12 12:22:21 1993
 *
 * CVS $Id: sqio.c,v 1.29 2002/08/26 23:10:52 eddy Exp)
 * 
 *****************************************************************
 * Basic API for single sequence reading:
 * 
 * SQFILE *sqfp;
 * char   *seqfile;
 * int     format;        - see squid.h for formats; example: SQFILE_FASTA
 * char   *seq;
 * SQINFO  sqinfo;
 * 
 * if ((sqfp = SeqfileOpen(seqfile, format, "BLASTDB")) == NULL)
 *   Die("Failed to open sequence database file %s\n%s\n", seqfile, usage);
 * while (ReadSeq(sqfp, sqfp->format, &seq, &sqinfo)) {
 *   do_stuff;
 *   FreeSequence(seq, &sqinfo);
 * }
 * SeqfileClose(sqfp);  
 * 
 *****************************************************************  
 */

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <ctype.h>

#ifndef SEEK_SET
#include <unistd.h>	
#endif

#include "squid.h"
#include "msa.h"
#include "ssi.h"

static void SeqfileGetLine(SQFILE *V);

#define kStartLength  500

static char *aminos      = "ABCDEFGHIKLMNPQRSTVWXYZ*";
static char *primenuc    = "ACGTUN";
static char *protonly    = "EFIPQZ";

static SQFILE *seqfile_open(char *filename, int format, char *env, int ssimode);

/* Function: SeqfileOpen()
 * 
 * Purpose : Open a sequence database file and prepare for reading
 *           sequentially.
 *           
 * Args:     filename - name of file to open
 *           format   - format of file
 *           env      - environment variable for path (e.g. BLASTDB)   
 *           ssimode  - -1, SSI_OFFSET_I32, or SSI_OFFSET_I64
 *
 *           Returns opened SQFILE ptr, or NULL on failure.
 */
SQFILE *
SeqfileOpen(char *filename, int format, char *env)
{
  return seqfile_open(filename, format, env, -1);
}
SQFILE *
SeqfileOpenForIndexing(char *filename, int format, char *env, int ssimode)
{
  return seqfile_open(filename, format, env, ssimode);
}
static SQFILE *
seqfile_open(char *filename, int format, char *env, int ssimode)
{
  SQFILE *dbfp;

  dbfp = (SQFILE *) MallocOrDie (sizeof(SQFILE));

  dbfp->ssimode  = ssimode;
  dbfp->rpl      = -1;		/* flag meaning "unset" */
  dbfp->lastrpl  = 0;
  dbfp->maxrpl   = 0;
  dbfp->bpl      = -1;		/* flag meaning "unset" */
  dbfp->lastbpl  = 0;
  dbfp->maxbpl   = 0;

  /* Open our file handle.
   * Three possibilities:
   *    1. normal file open
   *    2. filename = "-";    read from stdin
   *    3. filename = "*.gz"; read thru pipe from gzip 
   * If we're reading from stdin or a pipe, we can't reliably
   * back up, so we can't do two-pass parsers like the interleaved alignment   
   * formats.
   */
  if (strcmp(filename, "-") == 0)
    {
      dbfp->f         = stdin;
      dbfp->do_stdin  = TRUE; 
      dbfp->do_gzip   = FALSE;
      dbfp->fname     = sre_strdup("[STDIN]", -1);
    }
#ifndef SRE_STRICT_ANSI
  /* popen(), pclose() aren't portable to non-POSIX systems; disable */
  else if (Strparse("^.*\\.gz$", filename, 0))
    {
      char cmd[256];

      /* Note that popen() will return "successfully"
       * if file doesn't exist, because gzip works fine
       * and prints an error! So we have to check for
       * existence of file ourself.
       */
      if (! FileExists(filename))
	Die("%s: file does not exist", filename);

      if (strlen(filename) + strlen("gzip -dc ") >= 256)
	Die("filename > 255 char in SeqfileOpen()"); 
      sprintf(cmd, "gzip -dc %s", filename);
      if ((dbfp->f = popen(cmd, "r")) == NULL)
	return NULL;

      dbfp->do_stdin = FALSE;
      dbfp->do_gzip  = TRUE;
      dbfp->fname    = sre_strdup(filename, -1);
    }
#endif /*SRE_STRICT_ANSI*/
  else
    {
      if ((dbfp->f = fopen(filename, "r")) == NULL &&
	  (dbfp->f = EnvFileOpen(filename, env, NULL)) == NULL)
	return NULL;

      dbfp->do_stdin = FALSE;
      dbfp->do_gzip  = FALSE;
      dbfp->fname    = sre_strdup(filename, -1);
    }
  

  /* Invoke autodetection if we haven't already been told what
   * to expect.
   */
  if (format == SQFILE_UNKNOWN)
    {
      if (dbfp->do_stdin == TRUE || dbfp->do_gzip)
	Die("Can't autodetect sequence file format from a stdin or gzip pipe");
      format = SeqfileFormat(dbfp->f);
      if (format == SQFILE_UNKNOWN)
	Die("Can't determine format of sequence file %s", dbfp->fname);
    }

  /* The hack for sequential access of an interleaved alignment file:
   * read the alignment in, we'll copy sequences out one at a time.
   */
  dbfp->msa        = NULL;
  dbfp->afp        = NULL;
  dbfp->format     = format;
  dbfp->linenumber = 0;
  dbfp->buf        = NULL;
  dbfp->buflen     = 0;
  if (IsAlignmentFormat(format))	
    {
      /* We'll be reading from the MSA interface. Copy our data
       * to the MSA afp's structure.
       */
      dbfp->afp           = MallocOrDie(sizeof(MSAFILE));
      dbfp->afp->f        = dbfp->f;            /* just a ptr, don't close */
      dbfp->afp->do_stdin = dbfp->do_stdin;
      dbfp->afp->do_gzip  = dbfp->do_gzip;
      dbfp->afp->fname    = dbfp->fname;        /* just a ptr, don't free */
      dbfp->afp->format   = dbfp->format;       /* e.g. format */
      dbfp->afp->linenumber = dbfp->linenumber;	/* e.g. 0 */
      dbfp->afp->buf      = NULL;
      dbfp->afp->buflen   = 0;

      if ((dbfp->msa = MSAFileRead(dbfp->afp)) == NULL)
	Die("Failed to read any alignment data from file %s", dbfp->fname);
				/* hack: overload/reuse msa->lastidx; indicates
				   next seq to return upon a ReadSeq() call */
      dbfp->msa->lastidx = 0;

      return dbfp;
    }

  /* Load the first line.
   */
  SeqfileGetLine(dbfp); 
  return dbfp;
}

/* Function: SeqfilePosition()
 * 
 * Purpose:  Move to a particular offset in a seqfile.
 *           Will not work on alignment files.
 */
void
SeqfilePosition(SQFILE *sqfp, SSIOFFSET *offset)
{
  if (sqfp->do_stdin || sqfp->do_gzip || IsAlignmentFormat(sqfp->format))
    Die("SeqfilePosition() failed: in a nonrewindable data file or stream");

  if (SSISetFilePosition(sqfp->f, offset) != 0)
    Die("SSISetFilePosition failed, but that shouldn't happen.");
  SeqfileGetLine(sqfp);
}


/* Function: SeqfileRewind()
 * 
 * Purpose:  Set a sequence file back to the first sequence.
 * 
 *           Won't work on alignment files. Although it would
 *           seem that it could (just set msa->lastidx back to 0),
 *           that'll fail on "multiple multiple" alignment file formats
 *           (e.g. Stockholm).
 */
void
SeqfileRewind(SQFILE *sqfp)
{
  if (sqfp->do_stdin || sqfp->do_gzip)
    Die("SeqfileRewind() failed: in a nonrewindable data file or stream");

  rewind(sqfp->f);
  SeqfileGetLine(sqfp);
}

/* Function: SeqfileLineParameters()
 * Date:     SRE, Thu Feb 15 17:00:41 2001 [St. Louis]
 *
 * Purpose:  After all the sequences have been read from the file,
 *           but before closing it, retrieve overall bytes-per-line and
 *           residues-per-line info. If non-zero, these mean that
 *           the file contains homogeneous sequence line lengths (except
 *           the last line in each record).  
 *           
 *           If either of bpl or rpl is determined to be inhomogeneous,
 *           both are returned as 0.
 *
 * Args:     *sqfp   - an open but fully read sequence file
 *           ret_bpl - RETURN: bytes per line, or 0 if inhomogeneous
 *           ret_rpl - RETURN: residues per line, or 0 if inhomogenous.
 *
 * Returns:  void
 */
void
SeqfileLineParameters(SQFILE *V, int *ret_bpl, int *ret_rpl)
{
  if (V->rpl > 0 && V->maxrpl == V->rpl &&
      V->bpl > 0 && V->maxbpl == V->bpl) {
    *ret_bpl = V->bpl;
    *ret_rpl = V->rpl;
  } else {
    *ret_bpl = 0;
    *ret_rpl = 0;
  }
}


void
SeqfileClose(SQFILE *sqfp)
{
  /* note: don't test for sqfp->msa being NULL. Now that
   * we're holding afp open and allowing access to multi-MSA
   * databases (e.g. Stockholm format, Pfam), msa ends
   * up being NULL when we run out of alignments.
   */
  if (sqfp->afp != NULL) {
    if (sqfp->msa      != NULL) MSAFree(sqfp->msa);
    if (sqfp->afp->buf != NULL) free(sqfp->afp->buf);
    free(sqfp->afp);
  }
#ifndef SRE_STRICT_ANSI	/* gunzip functionality only on POSIX systems */
  if (sqfp->do_gzip)         pclose(sqfp->f);
#endif  
  else if (! sqfp->do_stdin) fclose(sqfp->f);
  if (sqfp->buf   != NULL) free(sqfp->buf);
  if (sqfp->fname != NULL) free(sqfp->fname);
  free(sqfp);
}


/* Function: SeqfileGetLine()
 * Date:     SRE, Tue Jun 22 09:15:49 1999 [Sanger Centre]
 *
 * Purpose:  read a line from a sequence file into V->buf
 *           If the fgets() is NULL, sets V->buf[0] to '\0'.
 *
 * Args:     V
 *
 * Returns:  void
 */
static void 
SeqfileGetLine(SQFILE *V)
{
  if (V->ssimode >= 0) 
    if (0 != SSIGetFilePosition(V->f, V->ssimode, &(V->ssioffset)))
      Die("SSIGetFilePosition() failed");
  if (sre_fgets(&(V->buf), &(V->buflen), V->f) == NULL)
    *(V->buf) = '\0';
  V->linenumber++;
}


void
FreeSequence(char *seq, SQINFO *sqinfo)
{
  if (seq != NULL) free(seq); /* FS, r244, here is potential problem in profile/profile */
  if (sqinfo->flags & SQINFO_SS){
    if (NULL != sqinfo->ss){ /* FS, r244 -> r245 */
      free(sqinfo->ss);
    }
  }
  if (sqinfo->flags & SQINFO_SA){
    if (NULL != sqinfo->sa){ /* FS, r244 -> r245 */
      free(sqinfo->sa);
    }
  }
}

int
SetSeqinfoString(SQINFO *sqinfo, char *sptr, int flag)
{
  int len;
  int pos;

				/* silently ignore NULL. */
  if (sptr == NULL) return 1;

  while (*sptr == ' ') sptr++; /* ignore leading whitespace */
  for (pos = strlen(sptr)-1; pos >= 0; pos--)
    if (! isspace((int) sptr[pos])) break;
  sptr[pos+1] = '\0';	       /* ignore trailing whitespace */

  switch (flag) {
  case SQINFO_NAME:
    if (*sptr != '-')
      { 
	strncpy(sqinfo->name, sptr, SQINFO_NAMELEN-1);
	sqinfo->name[SQINFO_NAMELEN-1] = '\0';
	sqinfo->flags   |= SQINFO_NAME;
      }
    break;

  case SQINFO_ID:
    if (*sptr != '-')
      { 
	strncpy(sqinfo->id, sptr, SQINFO_NAMELEN-1);
	sqinfo->id[SQINFO_NAMELEN-1] = '\0';
	sqinfo->flags |= SQINFO_ID;
      }
    break;

  case SQINFO_ACC:
    if (*sptr != '-')
      { 
	strncpy(sqinfo->acc, sptr, SQINFO_NAMELEN-1);
	sqinfo->acc[SQINFO_NAMELEN-1] = '\0';
	sqinfo->flags   |= SQINFO_ACC;
      }
    break;

  case SQINFO_DESC:
    if (*sptr != '-')
      { 
	if (sqinfo->flags & SQINFO_DESC) /* append? */
	  {
	    len = strlen(sqinfo->desc);
	    if (len < SQINFO_DESCLEN-2)	/* is there room? */
	      {
		strncat(sqinfo->desc, " ", SQINFO_DESCLEN-1-len); len++;
		strncat(sqinfo->desc, sptr, SQINFO_DESCLEN-1-len);
	      }
	  }
	else			/* else copy */
	  strncpy(sqinfo->desc, sptr, SQINFO_DESCLEN-1);
	sqinfo->desc[SQINFO_DESCLEN-1] = '\0';
	sqinfo->flags   |= SQINFO_DESC;
      }
    break;

  case SQINFO_START:
    if (!IsInt(sptr)) { squid_errno = SQERR_FORMAT; return 0; }
    sqinfo->start = atoi(sptr);
    if (sqinfo->start != 0) sqinfo->flags |= SQINFO_START;
    break;

  case SQINFO_STOP:
    if (!IsInt(sptr)) { squid_errno = SQERR_FORMAT; return 0; }
    sqinfo->stop = atoi(sptr);
    if (sqinfo->stop != 0) sqinfo->flags |= SQINFO_STOP;
    break;

  case SQINFO_OLEN:
    if (!IsInt(sptr)) { squid_errno = SQERR_FORMAT; return 0; }
    sqinfo->olen = atoi(sptr);
    if (sqinfo->olen != 0) sqinfo->flags |= SQINFO_OLEN;
    break;

  default:
    Die("Invalid flag %d to SetSeqinfoString()", flag);
  }
  return 1;
}

void
SeqinfoCopy(SQINFO *sq1, SQINFO *sq2)
{
  sq1->flags = sq2->flags;
  if (sq2->flags & SQINFO_NAME)  strcpy(sq1->name, sq2->name);
  if (sq2->flags & SQINFO_ID)    strcpy(sq1->id,   sq2->id);
  if (sq2->flags & SQINFO_ACC)   strcpy(sq1->acc,  sq2->acc);
  if (sq2->flags & SQINFO_DESC)  strcpy(sq1->desc, sq2->desc);
  if (sq2->flags & SQINFO_LEN)   sq1->len    = sq2->len;
  if (sq2->flags & SQINFO_START) sq1->start  = sq2->start;
  if (sq2->flags & SQINFO_STOP)  sq1->stop   = sq2->stop;
  if (sq2->flags & SQINFO_OLEN)  sq1->olen   = sq2->olen;
  if (sq2->flags & SQINFO_TYPE)  sq1->type   = sq2->type;
  if (sq2->flags & SQINFO_SS)    sq1->ss     = Strdup(sq2->ss);
  if (sq2->flags & SQINFO_SA)    sq1->sa     = Strdup(sq2->sa);
}

/* Function: ToDNA()
 * 
 * Purpose:  Convert a sequence to DNA.
 *           U --> T
 */
void
ToDNA(char *seq)
{
  for (; *seq != '\0'; seq++)
    {
      if      (*seq == 'U') *seq = 'T';
      else if (*seq == 'u') *seq = 't';
    }
}

/* Function: ToRNA()
 * 
 * Purpose:  Convert a sequence to RNA.
 *           T --> U
 */
void
ToRNA(char *seq)
{
  for (; *seq != '\0'; seq++)
    {
      if      (*seq == 'T') *seq = 'U';
      else if (*seq == 't') *seq = 'u';
    }
}


/* Function: ToIUPAC()
 * 
 * Purpose:  Convert X's, o's, other junk in a nucleic acid sequence to N's,
 *           to comply with IUPAC code. If is_aseq is TRUE, will allow gap
 *           characters though, so we can call ToIUPAC() on aligned seqs.
 *      
 *           NUCLEOTIDES is defined in squid.h as:
 *               "ACGTUNRYMKSWHBVDacgtunrymkswhbvd"
 *           gap chars allowed by isgap() are defined in squid.h as:
 *               " ._-~"
 *
 *           WU-BLAST's pressdb will
 *           choke on X's, for instance, necessitating conversion
 *           of certain genome centers' data. 
 */
void
ToIUPAC(char *seq, int is_aseq) 
{
  if (is_aseq) {
    for (; *seq != '\0'; seq++)
      if (strchr(NUCLEOTIDES, *seq) == NULL && ! isgap(*seq)) *seq = 'N';
  } else {
    for (; *seq != '\0'; seq++)
      if (strchr(NUCLEOTIDES, *seq) == NULL) *seq = 'N';
  }
}


/* Function: addseq()
 * 
 * Purpose:  Add a line of sequence to the growing string in V.
 *
 *           In the seven supported unaligned formats, all sequence
 *           lines may contain whitespace that must be filtered out;
 *           four formats (PIR, EMBL, Genbank, GCG) include coordinates
 *           that must be filtered out. Thus an (!isdigit && !isspace)
 *           test on each character before we accept it.
 */
static void 
addseq(char *s, struct ReadSeqVars *V)
{
  char *s0;
  char *sq;
  int   rpl;			/* valid residues per line */
  int   bpl;			/* characters per line     */

  if (V->ssimode == -1)
    {				/* Normal mode: keeping the seq */
      /* Make sure we have enough room. We know that s is <= buflen,
       * so just make sure we've got room for a whole new buflen worth
       * of sequence.
       */
      if (V->seqlen + V->buflen > V->maxseq) {
	V->maxseq += MAX(V->buflen, kStartLength);
	V->seq = ReallocOrDie (V->seq, V->maxseq+1);
      }

      sq = V->seq + V->seqlen;
      while (*s != 0) {
#ifdef CLUSTALO
    if (! isdigit((int) *s) && ! isspace((int) *s) && isprint((int) *s)) { 	
#else
	if (! isdigit((int) *s) && ! isspace((int) *s)) { 
#endif
	  *sq = *s;
	  sq++;
	}
	s++;
      }
      V->seqlen = sq - V->seq;
    }
  else				/* else: indexing mode, discard the seq */
    {
      s0 = s;
      rpl = 0;
      while (*s != 0) {
 	if (! isdigit((int) *s) && ! isspace((int) *s)) {
	  rpl++;
	}
	s++;
      }
      V->seqlen += rpl;
      bpl = s - s0;

      /* Keep track of the global rpl, bpl for the file.
       * This is overly complicated because we have to 
       * allow the last line of each record (e.g. the last addseq() call
       * on each sequence) to have a different length - and sometimes
       * we'll have one-line sequence records, too.  Thus we only
       * do something with the global V->rpl when we have *passed over*
       * a line - we keep the last line's rpl in last_rpl. And because
       * a file might consist entirely of single-line records, we keep
       * a third guy, maxrpl, that tells us the maximum rpl of any line
       * in the file. If we reach the end of file and rpl is still unset,
       * we'll set it to maxrpl. If we reach eof and rpl is set, but is
       * less than maxrpl, that's a weird case where a last line in some
       * record is longer than every other line.
       */
      if (V->rpl != 0) {		/* 0 means we already know rpl is invalid       */
	if (V->lastrpl > 0) {	/* we're on something that's not the first line */
	  if (V->rpl > 0 && V->lastrpl != V->rpl)  V->rpl = 0; 
	  else if (V->rpl == -1)                   V->rpl = V->lastrpl;
	}      
	V->lastrpl = rpl;
	if (rpl > V->maxrpl) V->maxrpl = rpl; /* make sure we check max length of final lines */
      }
      if (V->bpl != 0) {		/* 0 means we already know bpl is invalid       */
	if (V->lastbpl > 0) {	/* we're on something that's not the first line */
	  if (V->bpl > 0 && V->lastbpl != V->bpl)  V->bpl = 0; 
	  else if (V->bpl == -1)                   V->bpl = V->lastbpl;
	}      
	V->lastbpl = bpl;
	if (bpl > V->maxbpl) V->maxbpl = bpl; /* make sure we check max length of final lines */
      }
    } /* end of indexing mode of addseq(). */

}

static void 
readLoop(int addfirst, int (*endTest)(char *,int *), struct ReadSeqVars *V)
{
  int addend = 0;
  int done   = 0;

  V->seqlen = 0;
  V->lastrpl = V->lastbpl = 0;
  if (addfirst) {
    if (V->ssimode >= 0) V->d_off = V->ssioffset;
    addseq(V->buf, V);
  } else if (V->ssimode >= 0)
    if (0 != SSIGetFilePosition(V->f, V->ssimode, &(V->d_off)))
      Die("SSIGetFilePosition() failed");

  do {
    SeqfileGetLine(V);
	/* feof() alone is a bug; files not necessarily \n terminated */
    if (*(V->buf) == '\0' && feof(V->f))
      done = TRUE;
    done |= (*endTest)(V->buf, &addend);
    if (addend || !done)
      addseq(V->buf, V);
  } while (!done);
}


static int
endPIR(char *s, int  *addend)
{
  *addend = 0;
  if ((strncmp(s, "///", 3) == 0) || 
      (strncmp(s, "ENTRY", 5) == 0))
    return 1;
  else
    return 0;
}

static void
readPIR(struct ReadSeqVars *V)
{
  char *sptr;
				/* load first line of entry  */
  while (!feof(V->f) && strncmp(V->buf, "ENTRY", 5) != 0) {
    SeqfileGetLine(V);
  }
  if (feof(V->f)) return;
  if (V->ssimode >= 0) V->r_off = V->ssioffset;

  if ((sptr = strtok(V->buf + 15, "\n\t ")) != NULL)
    {
      SetSeqinfoString(V->sqinfo, sptr, SQINFO_NAME);
      SetSeqinfoString(V->sqinfo, sptr, SQINFO_ID);
    }
  do {
    SeqfileGetLine(V);
    if (!feof(V->f) && strncmp(V->buf, "TITLE", 5) == 0)
      SetSeqinfoString(V->sqinfo, V->buf+15, SQINFO_DESC);
    else if (!feof(V->f) && strncmp(V->buf, "ACCESSION", 9) == 0)
      {
	if ((sptr = strtok(V->buf+15, " \t\n")) != NULL)
	  SetSeqinfoString(V->sqinfo, sptr, SQINFO_ACC);
      }
  } while (! feof(V->f) && (strncmp(V->buf,"SEQUENCE", 8) != 0));
  SeqfileGetLine(V);			/* skip next line, coords */

  readLoop(0, endPIR, V);

  /* reading a real PIR-CODATA database file, we keep the source coords
   */
  V->sqinfo->start = 1;
  V->sqinfo->stop  = V->seqlen;
  V->sqinfo->olen  = V->seqlen;
  V->sqinfo->flags |= SQINFO_START | SQINFO_STOP | SQINFO_OLEN;

  /* get next line
   */
  while (!feof(V->f) && strncmp(V->buf, "ENTRY", 5) != 0) {
    SeqfileGetLine(V);
  }
}



static int 
endIG(char *s, int  *addend)
{
  *addend = 1; /* 1 or 2 occur in line w/ bases */
  return((strchr(s,'1')!=NULL) || (strchr(s,'2')!=NULL));
}

static void 
readIG(struct ReadSeqVars *V)
{
  char *nm;
				/* position past ';' comments */
  do {
    SeqfileGetLine(V);
  } while (! (feof(V->f) || ((*V->buf != 0) && (*V->buf != ';')) ));

  if (!feof(V->f))
    {
      if ((nm = strtok(V->buf, "\n\t ")) != NULL)
	SetSeqinfoString(V->sqinfo, nm, SQINFO_NAME);

      readLoop(0, endIG, V);
    }
  
  while (!(feof(V->f) || ((*V->buf != '\0') && (*V->buf == ';'))))
    SeqfileGetLine(V);
}

static int 
endStrider(char *s, int *addend)
{
  *addend = 0;
  return (strstr( s, "//") != NULL);
}

static void 
readStrider(struct ReadSeqVars *V)
{ 
  char *nm;
  
  while ((!feof(V->f)) && (*V->buf == ';')) 
    {
      if (strncmp(V->buf,"; DNA sequence", 14) == 0)
	{
	  if ((nm = strtok(V->buf+16, ",\n\t ")) != NULL)
	    SetSeqinfoString(V->sqinfo, nm, SQINFO_NAME);
	}
      SeqfileGetLine(V);
    }

  if (! feof(V->f))
    readLoop(1, endStrider, V);

  /* load next line
   */
  while ((!feof(V->f)) && (*V->buf != ';')) 
    SeqfileGetLine(V);
}


static int 
endGB(char *s, int *addend)
{
  *addend = 0;
  return ((strstr(s,"//") != NULL) || (strstr(s,"LOCUS") == s));
}

static void 
readGenBank(struct ReadSeqVars *V)
{
  char *sptr;
  int   in_definition;

  /* We'll map three genbank identifiers onto names:
   *     LOCUS     -> sqinfo.name
   *     ACCESSION -> sqinfo.acc   [primary accession only]
   *     VERSION   -> sqinfo.id
   * We don't currently store the GI number, or secondary accessions.    
   */
  while (strncmp(V->buf, "LOCUS", 5) != 0) {
    SeqfileGetLine(V);
  }
  if (V->ssimode >= 0) V->r_off = V->ssioffset;

  if ((sptr = strtok(V->buf+12, "\n\t ")) != NULL)
    SetSeqinfoString(V->sqinfo, sptr, SQINFO_NAME);

  in_definition = FALSE;
  while (! feof(V->f))
    {
      SeqfileGetLine(V);
      if (! feof(V->f) && strstr(V->buf, "DEFINITION") == V->buf)
	{
	  if ((sptr = strtok(V->buf+12, "\n")) != NULL)
	    SetSeqinfoString(V->sqinfo, sptr, SQINFO_DESC);
	  in_definition = TRUE;
	}
      else if (! feof(V->f) && strstr(V->buf, "ACCESSION") == V->buf)
	{
	  if ((sptr = strtok(V->buf+12, "\n\t ")) != NULL)
	    SetSeqinfoString(V->sqinfo, sptr, SQINFO_ACC);
	  in_definition = FALSE;
	}
      else if (! feof(V->f) && strstr(V->buf, "VERSION") == V->buf)
	{
	  if ((sptr = strtok(V->buf+12, "\n\t ")) != NULL)
	    SetSeqinfoString(V->sqinfo, sptr, SQINFO_ID);
	  in_definition = FALSE;
	}
      else if (strncmp(V->buf,"ORIGIN", 6) != 0)
	{
	  if (in_definition)
	    SetSeqinfoString(V->sqinfo, V->buf, SQINFO_DESC);
	}
      else
	break;
    }

  readLoop(0, endGB, V);

  /* reading a real GenBank database file, we keep the source coords
   */
  V->sqinfo->start = 1;
  V->sqinfo->stop  = V->seqlen;
  V->sqinfo->olen  = V->seqlen;
  V->sqinfo->flags |= SQINFO_START | SQINFO_STOP | SQINFO_OLEN;


  while (!(feof(V->f) || ((*V->buf!=0) && (strstr(V->buf,"LOCUS") == V->buf))))
    SeqfileGetLine(V);
				/* SRE: V->s now holds "//", so sequential
				   reads are wedged: fixed Tue Jul 13 1993 */
  while (!feof(V->f) && strstr(V->buf, "LOCUS  ") != V->buf)
    SeqfileGetLine(V);
}

static int
endGCGdata(char *s, int *addend) 
{
  *addend = 0;
  return (*s == '>');
}

static void
readGCGdata(struct ReadSeqVars *V)
{
  int   binary = FALSE;		/* whether data are binary or not */
  int   blen = 0;		/* length of binary sequence */
  
				/* first line contains ">>>>" followed by name */
  if (Strparse(">>>>([^ ]+) .+2BIT +Len: ([0-9]+)", V->buf, 2))
    {
      binary = TRUE;
      SetSeqinfoString(V->sqinfo, sqd_parse[1], SQINFO_NAME);
      blen = atoi(sqd_parse[2]);
    } 
  else if (Strparse(">>>>([^ ]+) .+ASCII +Len: [0-9]+", V->buf, 1))
    SetSeqinfoString(V->sqinfo, sqd_parse[1], SQINFO_NAME);
  else 
    Die("bogus GCGdata format? %s", V->buf);

				/* second line contains free text description */
  SeqfileGetLine(V);
  SetSeqinfoString(V->sqinfo, V->buf, SQINFO_DESC);

  if (binary) {
    /* allocate for blen characters +3... (allow for 3 bytes of slop) */
    if (blen >= V->maxseq) {
      V->maxseq = blen;
      if ((V->seq = (char *) realloc (V->seq, sizeof(char)*(V->maxseq+4)))==NULL)
	Die("malloc failed");
    }
				/* read (blen+3)/4 bytes from file */
    if (fread(V->seq, sizeof(char), (blen+3)/4, V->f) < (size_t) ((blen+3)/4))
      Die("fread failed");
    V->seqlen = blen;
				/* convert binary code to seq */
    GCGBinaryToSequence(V->seq, blen);
  }
  else readLoop(0, endGCGdata, V);
  
  while (!(feof(V->f) || ((*V->buf != 0) && (*V->buf == '>'))))
    SeqfileGetLine(V);
}

static int
endPearson(char *s, int *addend)
{
  *addend = 0;
  return(*s == '>');
}

static void 
readPearson(struct ReadSeqVars *V)
{
  char *sptr;

  if (V->ssimode >= 0) V->r_off = V->ssioffset;

  if (*V->buf != '>') 
    Die("\
File %s does not appear to be in FASTA format at line %d.\n\
You may want to specify the file format on the command line.\n\
Usually this is done with an option --informat <fmt>.\n", 
	V->fname, V->linenumber);

  if ((sptr = strtok(V->buf+1, "\n\t ")) != NULL)
    SetSeqinfoString(V->sqinfo, sptr, SQINFO_NAME);
  if ((sptr = strtok(NULL, "\n")) != NULL)
    SetSeqinfoString(V->sqinfo, sptr, SQINFO_DESC);

  readLoop(0, endPearson, V);

  while (!(feof(V->f) || ((*V->buf != 0) && (*V->buf == '>')))) {
    SeqfileGetLine(V);
  }
}


static int
endEMBL(char *s, int *addend)
{
  *addend = 0;
  /* Some people (Berlin 5S rRNA database, f'r instance) use
   * an extended EMBL format that attaches extra data after
   * the sequence -- watch out for that. We use the fact that
   * real EMBL sequence lines begin with five spaces.
   * 
   * We can use this as the sole end test because readEMBL() will
   * advance to the next ID line before starting to read again.
   */
  return (strncmp(s,"     ",5) != 0);
/*  return ((strstr(s,"//") != NULL) || (strstr(s,"ID   ") == s)); */
}

static void 
readEMBL(struct ReadSeqVars *V)
{
  char *sptr;

				/* make sure we have first line */
  while (!feof(V->f) && strncmp(V->buf, "ID  ", 4) != 0) {
    SeqfileGetLine(V);
  }
  if (V->ssimode >= 0) V->r_off = V->ssioffset;

  if ((sptr = strtok(V->buf+5, "\n\t ")) != NULL)
    {
      SetSeqinfoString(V->sqinfo, sptr, SQINFO_NAME);
      SetSeqinfoString(V->sqinfo, sptr, SQINFO_ID);
    }

  do {
    SeqfileGetLine(V);
    if (!feof(V->f) && strstr(V->buf, "AC  ") == V->buf)
      {
	if ((sptr = strtok(V->buf+5, ";  \t\n")) != NULL)
	  SetSeqinfoString(V->sqinfo, sptr, SQINFO_ACC);
      }
    else if (!feof(V->f) && strstr(V->buf, "DE  ") == V->buf)
      {
	if ((sptr = strtok(V->buf+5, "\n")) != NULL)
	  SetSeqinfoString(V->sqinfo, sptr, SQINFO_DESC);
      }
  } while (! feof(V->f) && strncmp(V->buf,"SQ",2) != 0);
  
  readLoop(0, endEMBL, V);

  /* Hack for Staden experiment files: convert - to N
   */
  if (V->ssimode == -1)		/* if we're in ssi mode, we're not keeping the seq */
    for (sptr = V->seq; *sptr != '\0'; sptr++)
      if (*sptr == '-') *sptr = 'N';

  /* reading a real EMBL database file, we keep the source coords
   */
  V->sqinfo->start = 1;
  V->sqinfo->stop  = V->seqlen;
  V->sqinfo->olen  = V->seqlen;
  V->sqinfo->flags |= SQINFO_START | SQINFO_STOP | SQINFO_OLEN;

				/* load next record's ID line */
  while (!feof(V->f) && strncmp(V->buf, "ID  ", 4) != 0) {
    SeqfileGetLine(V);
  }    

}


static int
endZuker(char *s, int *addend)
{
  *addend = 0;
  return( *s == '(' );
}

static void
readZuker(struct ReadSeqVars *V)
{
  char *sptr;

  SeqfileGetLine(V);  /*s == "seqLen seqid string..."*/

  if ((sptr = strtok(V->buf+6, " \t\n")) != NULL)
    SetSeqinfoString(V->sqinfo, sptr, SQINFO_NAME);

  if ((sptr = strtok(NULL, "\n")) != NULL)
    SetSeqinfoString(V->sqinfo, sptr, SQINFO_DESC);

  readLoop(0, endZuker, V);

  while (!(feof(V->f) | ((*V->buf != '\0') & (*V->buf == '('))))
    SeqfileGetLine(V);
}

static void 
readUWGCG(struct ReadSeqVars *V)
{
  char  *si;
  char  *sptr;
  int    done;

  V->seqlen = 0;

  /*writeseq: "    %s  Length: %d  (today)  Check: %d  ..\n" */
  /*drop above or ".." from id*/
  if ((si = strstr(V->buf,"  Length: ")) != NULL) *si = 0;
  else if ((si = strstr(V->buf,"..")) != NULL)    *si = 0;

  if ((sptr = strtok(V->buf, "\n\t ")) != NULL)
    SetSeqinfoString(V->sqinfo, sptr, SQINFO_NAME);

  do {
    done = feof(V->f);
    SeqfileGetLine(V);
    if (! done) addseq(V->buf, V);
  } while (!done);
}

    
/* Function: ReadSeq()
 * 
 * Purpose:  Read next sequence from an open database file.
 *           Return the sequence and associated info.
 *           
 * Args:     fp      - open sequence database file pointer          
 *           format  - format of the file (previously determined
 *                      by call to SeqfileFormat()).
 *                      Currently unused, since we carry it in V. 
 *           ret_seq - RETURN: sequence
 *           sqinfo  - RETURN: filled in w/ other information  
 *                     
 * Limitations: uses squid_errno, so it's not threadsafe.                    
 *           
 * Return:   1 on success, 0 on failure.
 *           ret_seq and some field of sqinfo are allocated here,
 *           The preferred call mechanism to properly free the memory is:
 *           
 *           SQINFO sqinfo;
 *           char  *seq;
 *           
 *           ReadSeq(fp, format, &seq, &sqinfo);
 *           ... do something...
 *           FreeSequence(seq, &sqinfo);
 */
int
ReadSeq(SQFILE *V, int format, char **ret_seq, SQINFO *sqinfo)
{
  int    gotuw;

  squid_errno = SQERR_OK;

  /* Here's the hack for sequential access of sequences from
   * the multiple sequence alignment formats
   */
  if (IsAlignmentFormat(V->format))
    {
      if (V->msa->lastidx >= V->msa->nseq) 
	{ /* out of data. try to read another alignment */				
	  MSAFree(V->msa);
	  if ((V->msa = MSAFileRead(V->afp)) == NULL)
	    return 0;
	  V->msa->lastidx = 0;
	}
				/* copy and dealign the appropriate aligned seq */
/* AW: stopping squid from dealigning sequences and corresponding info */
#ifdef CLUSTALO
      V->seq = sre_strdup(V->msa->aseq[V->msa->lastidx], V->msa->alen);
#else
      MakeDealignedString(V->msa->aseq[V->msa->lastidx], V->msa->alen, 
                          V->msa->aseq[V->msa->lastidx], &(V->seq));
#endif
      V->seqlen = strlen(V->seq);

      /* Extract sqinfo stuff for this sequence from the msa.
       * Tedious; code that should be cleaned.
       */
      sqinfo->flags = 0;
      if (V->msa->sqname[V->msa->lastidx] != NULL) 
	SetSeqinfoString(sqinfo, V->msa->sqname[V->msa->lastidx], SQINFO_NAME);
      if (V->msa->sqacc != NULL && V->msa->sqacc[V->msa->lastidx] != NULL) 
	SetSeqinfoString(sqinfo, V->msa->sqacc[V->msa->lastidx], SQINFO_ACC);
      if (V->msa->sqdesc != NULL && V->msa->sqdesc[V->msa->lastidx] != NULL) 
	SetSeqinfoString(sqinfo, V->msa->sqdesc[V->msa->lastidx], SQINFO_DESC);
      if (V->msa->ss != NULL && V->msa->ss[V->msa->lastidx] != NULL) {
/* AW: stopping squid from dealigning sequences and corresponding info */
#ifdef CLUSTALO
          sqinfo->ss = sre_strdup(V->msa->ss[V->msa->lastidx], V->msa->alen);
#else
          MakeDealignedString(V->msa->aseq[V->msa->lastidx], V->msa->alen, 
                              V->msa->ss[V->msa->lastidx], &(sqinfo->ss));
#endif
              sqinfo->flags |= SQINFO_SS;
      }
      if (V->msa->sa != NULL && V->msa->sa[V->msa->lastidx] != NULL) {
/* AW: stopping squid from dealigning sequences and corresponding info */
#ifdef CLUSTALO
          sqinfo->sa = sre_strdup(V->msa->sa[V->msa->lastidx], V->msa->alen);
#else
          MakeDealignedString(V->msa->aseq[V->msa->lastidx], V->msa->alen, 
                              V->msa->sa[V->msa->lastidx], &(sqinfo->sa));          
#endif
	sqinfo->flags |= SQINFO_SA;
      }
      V->msa->lastidx++;
    } 
  else {
    if (feof(V->f)) return 0;

    if (V->ssimode == -1) {	/* normal mode */
      V->seq           = (char*) calloc (kStartLength+1, sizeof(char));
      V->maxseq        = kStartLength;
    } else {			/* index mode: discarding seq */
      V->seq           = NULL;
      V->maxseq        = 0;
    }
    V->seqlen        = 0;
    V->sqinfo        = sqinfo;
    V->sqinfo->flags = 0;

    switch (V->format) {
    case SQFILE_IG      : readIG(V);      break;
    case SQFILE_STRIDER : readStrider(V); break;
    case SQFILE_GENBANK : readGenBank(V); break;
    case SQFILE_FASTA   : readPearson(V); break;
    case SQFILE_EMBL    : readEMBL(V);    break;
    case SQFILE_ZUKER   : readZuker(V);   break;
    case SQFILE_PIR     : readPIR(V);     break;
    case SQFILE_GCGDATA : readGCGdata(V); break; 
	
    case SQFILE_GCG :
      do {			/* skip leading comments on GCG file */
	gotuw = (strstr(V->buf,"..") != NULL);
	if (gotuw) readUWGCG(V);
	SeqfileGetLine(V);
      } while (! feof(V->f));
      break;

    case SQFILE_IDRAW:   /* SRE: no attempt to read idraw postscript */
    default:
      squid_errno = SQERR_FORMAT;
      free(V->seq);
      return 0;
    }
    if (V->seq != NULL)		/* (it can be NULL in indexing mode) */
      V->seq[V->seqlen] = 0; /* stick a string terminator on it */
  }

  /* Cleanup
   */
  sqinfo->len    = V->seqlen; 
  sqinfo->flags |= SQINFO_LEN;
  *ret_seq = V->seq;
  if (squid_errno == SQERR_OK) return 1; else return 0;
}  

/* Function: SeqfileFormat()
 * Date:     SRE, Tue Jun 22 10:58:58 1999 [Sanger Centre]
 *
 * Purpose:  Determine format of an open file.
 *           Returns format code.
 *           Rewinds the file.
 *           
 *           Autodetects the following unaligned formats:
 *              SQFILE_FASTA
 *              SQFILE_GENBANK
 *              SQFILE_EMBL
 *              SQFILE_GCG
 *              SQFILE_GCGDATA
 *              SQFILE_PIR
 *           Also autodetects the following alignment formats:
 *              MSAFILE_STOCKHOLM
 *              MSAFILE_MSF
 *              MSAFILE_CLUSTAL
 *              MSAFILE_SELEX
 *              MSAFILE_PHYLIP
 *
 *           Can't autodetect MSAFILE_A2M, calls it SQFILE_FASTA.
 *           MSAFileFormat() does the opposite.
 *
 * Args:     sfp -  open SQFILE
 *           
 * Return:   format code, or SQFILE_UNKNOWN if unrecognized
 */          
int
SeqfileFormat(FILE *fp)
{
  char *buf;
  int   len;
  int   fmt = SQFILE_UNKNOWN;
  int   ndataline;
  char *bufcpy, *s, *s1, *s2;
  int   has_junk;

  buf       = NULL;
  len       = 0;
  ndataline = 0;
  has_junk  = FALSE;
  while (sre_fgets(&buf, &len, fp) != NULL)
    {
      if (IsBlankline(buf)) continue;

      /* Well-behaved formats identify themselves in first nonblank line.
       */
      if (ndataline == 0)
	{
	  if (strncmp(buf, ">>>>", 4) == 0 && strstr(buf, "Len: "))
	    { fmt = SQFILE_GCGDATA; goto DONE; }

	  if (buf[0] == '>')
	    { fmt = SQFILE_FASTA; goto DONE; }

	  if (strncmp(buf, "!!AA_SEQUENCE", 13) == 0 ||
	      strncmp(buf, "!!NA_SEQUENCE", 13) == 0)
	    { fmt = SQFILE_GCG; goto DONE; }

	  if (strncmp(buf, "# STOCKHOLM 1.", 14) == 0)
	    { fmt = MSAFILE_STOCKHOLM; goto DONE; }

	  if (strncmp(buf, "CLUSTAL", 7) == 0 && 
	      strstr(buf, "multiple sequence alignment") != NULL)
	    { fmt = MSAFILE_CLUSTAL; goto DONE; }

	  if (strncmp(buf, "!!AA_MULTIPLE_ALIGNMENT", 23) == 0 ||
	      strncmp(buf, "!!NA_MULTIPLE_ALIGNMENT", 23) == 0)
	    { fmt = MSAFILE_MSF; goto DONE; }

	  			/* PHYLIP id: also just a good bet */
	  bufcpy = sre_strdup(buf, -1);
	  s = bufcpy;
	  if ((s1 = sre_strtok(&s, WHITESPACE, NULL)) != NULL &&
	      (s2 = sre_strtok(&s, WHITESPACE, NULL)) != NULL &&
	      IsInt(s1) && 
	      IsInt(s2))
	    { free(bufcpy); fmt = MSAFILE_PHYLIP; goto DONE; }
	  free(bufcpy);
	}

      /* We trust that other formats identify themselves soon.
       */
     				/* dead giveaways for extended SELEX */
      if (strncmp(buf, "#=AU", 4) == 0 ||
          strncmp(buf, "#=ID", 4) == 0 ||
	  strncmp(buf, "#=AC", 4) == 0 ||
	  strncmp(buf, "#=DE", 4) == 0 ||
	  strncmp(buf, "#=GA", 4) == 0 ||
	  strncmp(buf, "#=TC", 4) == 0 ||
	  strncmp(buf, "#=NC", 4) == 0 ||
	  strncmp(buf, "#=SQ", 4) == 0 ||
	  strncmp(buf, "#=SS", 4) == 0 ||
	  strncmp(buf, "#=CS", 4) == 0 ||
	  strncmp(buf, "#=RF", 4) == 0)
	{ fmt = MSAFILE_SELEX; goto DONE; }
	
      if (strncmp(buf, "///", 3) == 0 || strncmp(buf, "ENTRY ", 6) == 0)
	{ fmt = SQFILE_PIR; goto DONE; }

				/* a ha, diagnostic of an (old) MSF file */
      if ((strstr(buf, "..")    != NULL) && 
	  (strstr(buf, "MSF:")  != NULL) &&
	  (strstr(buf, "Check:")!= NULL))
	{ fmt = MSAFILE_MSF; goto DONE; }

				/* unaligned GCG (must follow MSF test!) */
      if (strstr(buf, " Check: ") != NULL && strstr(buf, "..") != NULL)
	{ fmt = SQFILE_GCG; goto DONE; }

      if (strncmp(buf,"LOCUS ",6) == 0 || strncmp(buf,"ORIGIN ",6) == 0)
	{ fmt = SQFILE_GENBANK; goto DONE; }

      if (strncmp(buf,"ID   ",5) == 0 || strncmp(buf,"SQ   ",5) == 0)
	{ fmt = SQFILE_EMBL; goto DONE; }

      /* But past here, we're being desperate. A simple SELEX file is
       * very difficult to detect; we can only try to disprove it.
       */
      s = buf;
      if ((s1 = sre_strtok(&s, WHITESPACE, NULL)) == NULL) continue; /* skip blank lines */
      if (strchr("#%", *s1) != NULL) continue;   /* skip comment lines */

      /* Disproof 1. Noncomment, nonblank lines in a SELEX file
       * must have at least two space-delimited fields (name/seq)
       */
      if ((s2 = sre_strtok(&s, WHITESPACE, NULL)) == NULL) 
	has_junk = TRUE;

      /* Disproof 2. 
       * The sequence field should look like a sequence.
       */
      if (s2 != NULL && Seqtype(s2) == kOtherSeq) 
	has_junk = TRUE;

      ndataline++;
      if (ndataline == 300) break; /* only look at first 300 lines */
    }

  if (ndataline == 0)
    Die("Sequence file contains no data");

  /* If we've made it this far, we've run out of data, but there
   * was at least one line of it; check if we've
   * disproven SELEX. If not, cross our fingers, pray, and guess SELEX. 
   */
  if (has_junk == TRUE) fmt = SQFILE_UNKNOWN;
  else                  fmt = MSAFILE_SELEX;

 DONE:
  if (buf != NULL) free(buf);
  rewind(fp);
  return fmt;
}

/* Function: GCGBinaryToSequence()
 * 
 * Purpose:  Convert a GCG 2BIT binary string to DNA sequence.
 *           0 = C  1 = T  2 = A  3 = G
 *           4 nts/byte
 *           
 * Args:     seq  - binary sequence. Converted in place to DNA.
 *           len  - length of DNA. binary is (len+3)/4 bytes
 */
int
GCGBinaryToSequence(char *seq, int len)
{
  int   bpos;			/* position in binary   */
  int   spos;			/* position in sequence */
  char  twobit;
  int   i;

  for (bpos = (len-1)/4; bpos >= 0; bpos--) 
    {
      twobit = seq[bpos];
      spos   = bpos*4;

      for (i = 3; i >= 0; i--) 
	{
	  switch (twobit & 0x3) {
	  case 0: seq[spos+i] = 'C'; break;
	  case 1: seq[spos+i] = 'T'; break;
	  case 2: seq[spos+i] = 'A'; break;
	  case 3: seq[spos+i] = 'G'; break;
	  }
	  twobit = twobit >> 2;
	}
    }
  seq[len] = '\0';
  return 1;
}


/* Function: GCGchecksum()
 * Date:     SRE, Mon May 31 11:13:21 1999 [St. Louis]
 *
 * Purpose:  Calculate a GCG checksum for a sequence.
 *           Code provided by Steve Smith of Genetics
 *           Computer Group.
 *
 * Args:     seq -  sequence to calculate checksum for.
 *                  may contain gap symbols.
 *           len -  length of sequence (usually known,
 *                  so save a strlen() call)       
 *
 * Returns:  GCG checksum.
 */
int
GCGchecksum(char *seq, int len)
{
  int i;			/* position in sequence */
  int chk = 0;			/* calculated checksum  */

  for (i = 0; i < len; i++) 
    chk = (chk + (i % 57 + 1) * (sre_toupper((int) seq[i]))) % 10000;
  return chk;
}


/* Function: GCGMultchecksum()
 * 
 * Purpose:  GCG checksum for a multiple alignment: sum of
 *           individual sequence checksums (including their
 *           gap characters) modulo 10000.
 *
 *           Implemented using spec provided by Steve Smith of
 *           Genetics Computer Group.
 *           
 * Args:     seqs - sequences to be checksummed; aligned or not
 *           nseq - number of sequences
 *           
 * Return:   the checksum, a number between 0 and 9999
 */                      
int
GCGMultchecksum(char **seqs, int nseq)
{
  int chk = 0;
  int idx;

  for (idx = 0; idx < nseq; idx++)
    chk = (chk + GCGchecksum(seqs[idx], strlen(seqs[idx]))) % 10000;
  return chk;
}




/* Function: Seqtype()
 * 
 * Purpose:  Returns a (very good) guess about type of sequence:
 *           kDNA, kRNA, kAmino, or kOtherSeq.
 *           
 *           Modified from, and replaces, Gilbert getseqtype().
 */
int
Seqtype(char *seq)
{
  int  saw;			/* how many non-gap characters I saw */
  char c;
  int  po = 0;			/* count of protein-only */
  int  nt = 0;			/* count of t's */
  int  nu = 0;			/* count of u's */
  int  na = 0;			/* count of nucleotides */
  int  aa = 0;			/* count of amino acids */
  int  no = 0;			/* count of others */
  
  /* Look at the first 300 non-gap characters
   */

#ifdef CLUSTALO
  /* VGGNGDDYLSGGTGNDTL is recognized as unknown using squid's default
   * approach. 
   * We change it to the following:

   * 1. counting: ignore gaps and not alpha characters. if protein-only then
   * count as such (po). otherwise decide if amino-acid (aa) or nucleic-acid
   * (na) or unknown (no)
   *
   * 2. determine type: if we saw more unknown than aa or na, return unknown.
   * if encountered protein-only return protein-only. otherwise decide based
   * on majority. (if aa==na return na)
   */
  for (saw = 0; *seq != '\0' && saw < 300; seq++) {
		  c = sre_toupper((int) *seq);
		  int unknown = 1;

		  if (isgap(c) ||  ! isalpha((int) c))  {
				  continue;
		  }

		  if (strchr(protonly, c)) {
				  po++;
				  unknown = 0;
		  }

		  if (strchr(aminos,c)) {
				  aa++;
				  unknown = 0;
		  }

		  if (strchr(primenuc,c)) {
				  na++;
				  unknown = 0;

				  if (c == 'T') 
						  nt++;
				  else if (c == 'U') 
						  nu++;
		  } 

		  if (unknown) {
				  no ++;
		  }

		  saw++;
  }

  if (no > aa && no > na)
		  return kOtherSeq;

 if (po > 0 || aa>na) 
		  return kAmino;

 if (na >= aa) {
		 if (nu > nt) 
				 return kRNA;
		 else 
				 return kDNA;
  }

  return kOtherSeq;


#else
  for (saw = 0; *seq != '\0' && saw < 300; seq++)
    {
      c = sre_toupper((int) *seq);
      if (! isgap(c)) 
	{
	  if (strchr(protonly, c)) po++;
	  else if (strchr(primenuc,c)) {
	    na++;
	    if (c == 'T') nt++;
	    else if (c == 'U') nu++;
	  }
	  else if (strchr(aminos,c)) aa++;
	  else if (isalpha((int) c)) no++;
	  saw++;
	}
    }

  if (no > 0) return kOtherSeq;
  else if (po > 0) return kAmino;
  else if (na > aa) {
    if (nu > nt) return kRNA;
    else return kDNA;
    }
  else return kAmino;		/* ooooh. risky. */
#endif

}


/* Function: GuessAlignmentSeqtype()
 * Date:     SRE, Wed Jul  7 09:42:34 1999 [St. Louis]
 *
 * Purpose:  Try to guess whether an alignment is protein 
 *           or nucleic acid; return a code for the
 *           type (kRNA, kDNA, or kAmino).
 *
 * Args:     aseq  - array of aligned sequences. (Could also
 *                   be an rseq unaligned sequence array)
 *           nseq  - number of aseqs
 *
 * Returns:  kRNA, kDNA, kAmino;
 *           kOtherSeq if inconsistency is detected.
 */
int
GuessAlignmentSeqtype(char **aseq, int nseq)
{
  int idx;
  int nrna   = 0;
  int ndna   = 0;
  int namino = 0;
  int nother = 0;

  for (idx = 0; idx < nseq; idx++)
    switch (Seqtype(aseq[idx])) {
    case kRNA:   nrna++;   break;
    case kDNA:   ndna++;   break;
    case kAmino: namino++; break;
    default:     nother++;
    }

  /* Unambiguous decisions:
   */
  if (nother)         return kOtherSeq;
  if (namino == nseq) return kAmino;
  if (ndna   == nseq) return kDNA;
  if (nrna   == nseq) return kRNA;

  /* Ambiguous decisions:
   */
  if (namino == 0)    return kRNA; /* it's nucleic acid, but seems mixed RNA/DNA */
  return kAmino;		   /* some amino acid seen; others probably short seqs, some 
				      of which may be entirely ACGT (ala,cys,gly,thr). We 
				      could be a little more sophisticated: U would be a giveaway
				      that we're not in protein seqs */
}

/* Function: WriteSimpleFASTA()
 * Date:     SRE, Tue Nov 16 18:06:00 1999 [St. Louis]
 *
 * Purpose:  Just write a FASTA format sequence to a file;
 *           minimal interface, mostly for quick and dirty programs.
 *
 * Args:     fp   - open file handle (stdout, possibly)
 *           seq  - sequence to output
 *           name - name for the sequence
 *           desc - optional description line, or NULL.
 *
 * Returns:  void
 */
void
WriteSimpleFASTA(FILE *fp, char *seq, char *name, char *desc)
{
  char buf[61];
  int  len;
  int  pos;
  
  len = strlen(seq);
  buf[60] = '\0';
  fprintf(fp, ">%s %s\n", name, desc != NULL ? desc : "");
  for (pos = 0; pos < len; pos += 60)
    {
      strncpy(buf, seq+pos, 60);
      fprintf(fp, "%s\n", buf);
    }
}

int
WriteSeq(FILE *outf, int outform, char *seq, SQINFO *sqinfo)
{
  int   numline = 0;
  int   lines = 0, spacer = 0, width  = 50, tab = 0;
  int   i, j, l, l1, ibase;
  char  endstr[10];
  char  s[100];			/* buffer for sequence  */
  char  ss[100];		/* buffer for structure */
  int   checksum = 0;
  int   seqlen;   
  int   which_case;    /* 0 = do nothing. 1 = upper case. 2 = lower case */
  int   dostruc;		/* TRUE to print structure lines*/

  which_case = 0;
  dostruc    = FALSE;		
  seqlen     = (sqinfo->flags & SQINFO_LEN) ? sqinfo->len : strlen(seq);

  if (IsAlignmentFormat(outform)) 
    Die("Tried to write an aligned format with WriteSeq() -- bad, bad.");


  strcpy( endstr,"");
  l1 = 0;
  checksum = GCGchecksum(seq, seqlen);

  switch (outform) {
  case SQFILE_UNKNOWN:    /* no header, just sequence */
    strcpy(endstr,"\n"); /* end w/ extra blank line */
    break;

  case SQFILE_GENBANK:
    fprintf(outf,"LOCUS       %s       %d bp\n", 
	    sqinfo->name, seqlen);
    fprintf(outf,"ACCESSION   %s\n", 
	    (sqinfo->flags & SQINFO_ACC) ? sqinfo->acc : ".");
    fprintf(outf,"DEFINITION  %s\n", 
	    (sqinfo->flags & SQINFO_DESC) ? sqinfo->desc : ".");
    fprintf(outf,"VERSION     %s\n", 
	    (sqinfo->flags & SQINFO_ID) ? sqinfo->id : ".");
    fprintf(outf,"ORIGIN      \n");
    spacer = 11;
    numline = 1;
    strcpy(endstr, "\n//");
    break;

  case SQFILE_GCGDATA:
    fprintf(outf, ">>>>%s  9/95  ASCII  Len: %d\n", sqinfo->name, seqlen);
    fprintf(outf, "%s\n", (sqinfo->flags & SQINFO_DESC) ? sqinfo->desc : "-");
    break;

  case SQFILE_PIR:
    fprintf(outf, "ENTRY          %s\n", 
	    (sqinfo->flags & SQINFO_ID) ? sqinfo->id : sqinfo->name);
    fprintf(outf, "TITLE          %s\n", 
	    (sqinfo->flags & SQINFO_DESC) ? sqinfo->desc : "-");
    fprintf(outf, "ACCESSION      %s\n",
	    (sqinfo->flags & SQINFO_ACC) ? sqinfo->acc : "-");
    fprintf(outf, "SUMMARY                                #Length %d  #Checksum  %d\n",
	    sqinfo->len, checksum);
    fprintf(outf, "SEQUENCE\n");
    fprintf(outf, "                  5        10        15        20        25        30\n");
    spacer  = 2;		/* spaces after every residue */
    numline = 1;              /* number lines w/ coords     */
    width   = 30;             /* 30 aa per line             */
    strcpy(endstr, "\n///");
    break;

  case SQFILE_SQUID:
    fprintf(outf, "NAM  %s\n", sqinfo->name);
    if (sqinfo->flags & (SQINFO_ID | SQINFO_ACC | SQINFO_START | SQINFO_STOP | SQINFO_OLEN))
      fprintf(outf, "SRC  %s %s %d..%d::%d\n",
	      (sqinfo->flags & SQINFO_ID)    ? sqinfo->id     : "-",
	      (sqinfo->flags & SQINFO_ACC)   ? sqinfo->acc    : "-",
	      (sqinfo->flags & SQINFO_START) ? sqinfo->start  : 0,
	      (sqinfo->flags & SQINFO_STOP)  ? sqinfo->stop   : 0,
	      (sqinfo->flags & SQINFO_OLEN)  ? sqinfo->olen   : 0);
    if (sqinfo->flags & SQINFO_DESC)
      fprintf(outf, "DES  %s\n", sqinfo->desc);
    if (sqinfo->flags & SQINFO_SS)
      {
	fprintf(outf, "SEQ  +SS\n");
	dostruc = TRUE;	/* print structure lines too */
      }
    else
      fprintf(outf, "SEQ\n");
    numline = 1;                /* number seq lines w/ coords  */
    strcpy(endstr, "\n++");
    break;

  case SQFILE_EMBL:
    fprintf(outf,"ID   %s\n",
	    (sqinfo->flags & SQINFO_ID) ? sqinfo->id : sqinfo->name);
    fprintf(outf,"AC   %s\n",
	    (sqinfo->flags & SQINFO_ACC) ? sqinfo->acc : "-");
    fprintf(outf,"DE   %s\n", 
	    (sqinfo->flags & SQINFO_DESC) ? sqinfo->desc : "-");
    fprintf(outf,"SQ             %d BP\n", seqlen);
    strcpy(endstr, "\n//"); /* 11Oct90: bug fix*/
    tab = 5;     /** added 31jan91 */
    spacer = 11; /** added 31jan91 */
    break;

  case SQFILE_GCG:
    fprintf(outf,"%s\n", sqinfo->name);
    if (sqinfo->flags & SQINFO_ACC)
      fprintf(outf,"ACCESSION   %s\n", sqinfo->acc); 
    if (sqinfo->flags & SQINFO_DESC)
      fprintf(outf,"DEFINITION  %s\n", sqinfo->desc);
    fprintf(outf,"    %s  Length: %d  (today)  Check: %d  ..\n", 
	    sqinfo->name, seqlen, checksum);
    spacer = 11;
    numline = 1;
    strcpy(endstr, "\n");  /* this is insurance to help prevent misreads at eof */
    break;

  case SQFILE_STRIDER: /* ?? map ?*/
    fprintf(outf,"; ### from DNA Strider ;-)\n");
    fprintf(outf,"; DNA sequence  %s, %d bases, %d checksum.\n;\n", 
	    sqinfo->name, seqlen, checksum);
    strcpy(endstr, "\n//");
    break;

    /* SRE: Don had Zuker default to Pearson, which is not
			   intuitive or helpful, since Zuker's MFOLD can't read
			   Pearson format. More useful to use kIG */
  case SQFILE_ZUKER:
    which_case = 1;			/* MFOLD requires upper case. */
    /*FALLTHRU*/
  case SQFILE_IG:
    fprintf(outf,";%s %s\n", 
	    sqinfo->name,
	    (sqinfo->flags & SQINFO_DESC) ? sqinfo->desc : "");
    fprintf(outf,"%s\n", sqinfo->name);
    strcpy(endstr,"1"); /* == linear dna */
    break;

  case SQFILE_RAW:			/* Raw: no header at all. */
    break;

  default :
  case SQFILE_FASTA:
    fprintf(outf,">%s %s\n", sqinfo->name,
	    (sqinfo->flags & SQINFO_DESC)  ? sqinfo->desc   : "");
    break;
  }

  if (which_case == 1) s2upper(seq);
  if (which_case == 2) s2lower(seq);


  width = MIN(width,100);
  for (i=0, l=0, ibase = 1, lines = 0; i < seqlen; ) {
    if (l1 < 0) l1 = 0;
    else if (l1 == 0) {
      if (numline) fprintf(outf,"%8d ",ibase);
      for (j=0; j<tab; j++) fputc(' ',outf);
    }
    if ((spacer != 0) && ((l+1) % spacer == 1)) 
      { s[l] = ' '; ss[l] = ' '; l++; }
    s[l]  = seq[i];
    ss[l] = (sqinfo->flags & SQINFO_SS) ? sqinfo->ss[i] : '.';
    l++; i++;
    l1++;                 /* don't count spaces for width*/
    if (l1 == width || i == seqlen) {
      s[l] = ss[l] = '\0';
      l = 0; l1 = 0;
      if (dostruc)
	{
	  fprintf(outf, "%s\n", s);
	  if (numline) fprintf(outf,"         ");
	  for (j=0; j<tab; j++) fputc(' ',outf);
	  if (i == seqlen) fprintf(outf,"%s%s\n",ss,endstr);
	  else fprintf(outf,"%s\n",ss);
	}
      else
	{
	  if (i == seqlen) fprintf(outf,"%s%s\n",s,endstr);
	  else fprintf(outf,"%s\n",s);
	}
      lines++;
      ibase = i+1;
    }
  }
  return lines;
} 


/* Function: ReadMultipleRseqs()
 * 
 * Purpose:  Open a data file and
 *           parse it into an array of rseqs (raw, unaligned
 *           sequences).
 * 
 *           Caller is responsible for free'ing memory allocated
 *           to ret_rseqs, ret_weights, and ret_names.
 *           
 *           Weights are currently only supported for MSF format.
 *           Sequences read from all other formats will be assigned
 *           weights of 1.0. If the caller isn't interested in
 *           weights, it passes NULL as ret_weights.
 * 
 * Returns 1 on success. Returns 0 on failure and sets
 * squid_errno to indicate the cause.
 */
int
ReadMultipleRseqs(char              *seqfile,
		  int                fformat,
		  char            ***ret_rseqs,
		  SQINFO **ret_sqinfo,
		  int               *ret_num)
{
  SQINFO *sqinfo;               /* array of sequence optional info         */
  SQFILE *dbfp;                 /* open ptr for sequential access of file  */
  char  **rseqs;                /* sequence array                          */
  int     numalloced;           /* num of seqs currently alloced for       */
  int     num;


  num        = 0;
  numalloced = 16;
  rseqs  = (char **) MallocOrDie (numalloced * sizeof(char *));
  sqinfo = (SQINFO *) MallocOrDie (numalloced * sizeof(SQINFO));
  if ((dbfp = SeqfileOpen(seqfile, fformat, NULL)) == NULL) return 0;      

  while (ReadSeq(dbfp, dbfp->format, &rseqs[num], &(sqinfo[num])))
    {
      num++;
      if (num == numalloced) /* more seqs coming, alloc more room */
	{
	  numalloced += 16;
	  rseqs  = (char **) ReallocOrDie (rseqs, numalloced*sizeof(char *));
	  sqinfo = (SQINFO *) ReallocOrDie (sqinfo, numalloced * sizeof(SQINFO));
	}
    }
  SeqfileClose(dbfp);

  *ret_rseqs  = rseqs;
  *ret_sqinfo = sqinfo;
  *ret_num    = num;
  return 1;
}


/* Function: String2SeqfileFormat()
 * Date:     SRE, Sun Jun 27 15:25:54 1999 [TW 723 over Canadian Shield]
 *
 * Purpose:  Convert a string (e.g. from command line option arg)
 *           to a format code. Case insensitive. Return
 *           MSAFILE_UNKNOWN/SQFILE_UNKNOWN if string is bad.  
 *           Uses codes defined in squid.h (unaligned formats) and
 *           msa.h (aligned formats).
 *
 * Args:     s   - string to convert; e.g. "stockholm"
 *
 * Returns:  format code; e.g. MSAFILE_STOCKHOLM
 */
int
String2SeqfileFormat(char *s)
{
  char *s2;
  int   code = SQFILE_UNKNOWN;

  if (s == NULL) return SQFILE_UNKNOWN;
  s2 = sre_strdup(s, -1);
  s2upper(s2);
  
  if      (strcmp(s2, "FASTA")     == 0) code = SQFILE_FASTA;
#ifdef CLUSTALO
  if      (strcmp(s2, "FA")        == 0) code = SQFILE_FASTA;
  else if (strcmp(s2, "VIENNA")    == 0) code = SQFILE_VIENNA;
  else if (strcmp(s2, "VIE")       == 0) code = SQFILE_VIENNA;
#endif
  else if (strcmp(s2, "GENBANK")   == 0) code = SQFILE_GENBANK;
#ifdef CLUSTALO
  else if (strcmp(s2, "GB")   == 0) code = SQFILE_GENBANK;
#endif
  else if (strcmp(s2, "EMBL")      == 0) code = SQFILE_EMBL;
  else if (strcmp(s2, "GCG")       == 0) code = SQFILE_GCG;
  else if (strcmp(s2, "GCGDATA")   == 0) code = SQFILE_GCGDATA;
  else if (strcmp(s2, "RAW")       == 0) code = SQFILE_RAW;
  else if (strcmp(s2, "IG")        == 0) code = SQFILE_IG;
  else if (strcmp(s2, "STRIDER")   == 0) code = SQFILE_STRIDER;
  else if (strcmp(s2, "IDRAW")     == 0) code = SQFILE_IDRAW;
  else if (strcmp(s2, "ZUKER")     == 0) code = SQFILE_ZUKER;
  else if (strcmp(s2, "PIR")       == 0) code = SQFILE_PIR;
  else if (strcmp(s2, "SQUID")     == 0) code = SQFILE_SQUID;
  else if (strcmp(s2, "STOCKHOLM") == 0) code = MSAFILE_STOCKHOLM;
#ifdef CLUSTALO
  else if (strcmp(s2, "ST") == 0) code = MSAFILE_STOCKHOLM;
  else if (strcmp(s2, "STK") == 0) code = MSAFILE_STOCKHOLM;
#endif
  else if (strcmp(s2, "SELEX")     == 0) code = MSAFILE_SELEX; 
  else if (strcmp(s2, "MSF")       == 0) code = MSAFILE_MSF; 
  else if (strcmp(s2, "CLUSTAL")   == 0) code = MSAFILE_CLUSTAL; 
#ifdef CLUSTALO
  else if (strcmp(s2, "CLU")   == 0) code = MSAFILE_CLUSTAL; 
#endif
  else if (strcmp(s2, "A2M")       == 0) code = MSAFILE_A2M; 
  else if (strcmp(s2, "PHYLIP")    == 0) code = MSAFILE_PHYLIP; 
#ifdef CLUSTALO
  else if (strcmp(s2, "PHY")    == 0) code = MSAFILE_PHYLIP; 
#endif
  else if (strcmp(s2, "EPS")       == 0) code = MSAFILE_EPS; 
#ifdef CLUSTALO
  else code = SQFILE_UNKNOWN;
#endif
  free(s2);
  return code;
}
char *
SeqfileFormat2String(int code)
{
  switch (code) {
  case SQFILE_UNKNOWN:   return "unknown";
  case SQFILE_FASTA:     return "FASTA";
#ifdef CLUSTALO
  case SQFILE_VIENNA:    return "Vienna";
#endif
  case SQFILE_GENBANK:   return "Genbank";
  case SQFILE_EMBL:      return "EMBL"; 
  case SQFILE_GCG:       return "GCG";
  case SQFILE_GCGDATA:   return "GCG data library";
  case SQFILE_RAW:       return "raw"; 
  case SQFILE_IG:        return "Intelligenetics";
  case SQFILE_STRIDER:   return "MacStrider";
  case SQFILE_IDRAW:     return "Idraw Postscript";
  case SQFILE_ZUKER:     return "Zuker"; 
  case SQFILE_PIR:       return "PIR";
  case SQFILE_SQUID:     return "SQUID";
  case MSAFILE_STOCKHOLM: return "Stockholm";
  case MSAFILE_SELEX:     return "SELEX";
  case MSAFILE_MSF:       return "MSF";
  case MSAFILE_CLUSTAL:   return "Clustal";
  case MSAFILE_A2M:       return "a2m";
  case MSAFILE_PHYLIP:    return "Phylip";
  case MSAFILE_EPS:       return "EPS";
  default:               
    Die("Bad code passed to MSAFormat2String()");
  }
  /*NOTREACHED*/
  return NULL;
}


/* Function: MSAToSqinfo()
 * Date:     SRE, Tue Jul 20 14:36:56 1999 [St. Louis]
 *
 * Purpose:  Take an MSA and generate a SQINFO array suitable
 *           for use in annotating the unaligned sequences.
 *           Return the array.
 *           
 *           Permanent temporary code. sqinfo was poorly designed.
 *           it must eventually be replaced, but the odds
 *           of this happening soon are nil, so I have to deal.
 *
 * Args:     msa   - the alignment
 *
 * Returns:  ptr to allocated sqinfo array.
 *           Freeing is ghastly: free in each individual sqinfo[i] 
 *           with FreeSequence(NULL, &(sqinfo[i])), then
 *           free(sqinfo).
 */
SQINFO *
MSAToSqinfo(MSA *msa)
{
  int     idx;
  SQINFO *sqinfo;

  sqinfo = MallocOrDie(sizeof(SQINFO) * msa->nseq);

  for (idx = 0; idx < msa->nseq; idx++)
    {
      sqinfo[idx].flags = 0;
      SetSeqinfoString(&(sqinfo[idx]), 
		       msa->sqname[idx],                 SQINFO_NAME);
      SetSeqinfoString(&(sqinfo[idx]), 
		       MSAGetSeqAccession(msa, idx),     SQINFO_ACC);
      SetSeqinfoString(&(sqinfo[idx]), 
		       MSAGetSeqDescription(msa, idx),   SQINFO_DESC);

      if (msa->ss != NULL && msa->ss[idx] != NULL) {
	MakeDealignedString(msa->aseq[idx], msa->alen, 
			    msa->ss[idx], &(sqinfo[idx].ss));
	sqinfo[idx].flags |= SQINFO_SS;
      }

      if (msa->sa != NULL && msa->sa[idx] != NULL) {
	MakeDealignedString(msa->aseq[idx], msa->alen, 
			    msa->sa[idx], &(sqinfo[idx].sa));
	sqinfo[idx].flags |= SQINFO_SA;
      }

      sqinfo[idx].len    = DealignedLength(msa->aseq[idx]);
      sqinfo[idx].flags |= SQINFO_LEN;
    }
  return sqinfo;
}



/* cc -o sqio_test -DA_QUIET_DAY -L. sqio.c -lsquid */
#ifdef A_QUIET_DAY
#include "ssi.h"
int
main(int argc, char **argv)
{
  FILE *fp;
  char *filename;
  char *buf;
  int   len;
  int   mode = 3;
  SSIOFFSET off;

  filename = argv[1];

  if (mode == 1) {
    buf = malloc(sizeof(char) * 256);
    if ((fp = fopen(filename, "r")) == NULL)
      Die("open of %s failed", filename); 
    while (fgets(buf, 255, fp) != NULL)
      ;
    fclose(fp);
    free(buf);
  } else if (mode == 2) {
    if ((fp = fopen(filename, "r")) == NULL)
      Die("open of %s failed", filename); 
    buf = NULL; len = 0;
    while (sre_fgets(&buf, &len, fp) != NULL)
      SSIGetFilePosition(fp, SSI_OFFSET_I32, &off); 
    fclose(fp);
    free(buf);
  } else if (mode == 3) {
    SQFILE *dbfp;
    SQINFO  info;

    if ((dbfp = SeqfileOpen(filename, SQFILE_FASTA, NULL)) == NULL)
      Die("open of %s failed", filename); 
    while (ReadSeq(dbfp, dbfp->format, &buf, &info)) { 
      SSIGetFilePosition(dbfp->f, SSI_OFFSET_I32, &off); 
      FreeSequence(buf, &info);
    }
    SeqfileClose(dbfp);
  }
      
}
 

#endif