/* seq2mtx - convert single sequence to pseudo IMPALA mtx file */

/* Copyright (C) 2000 D.T. Jones */

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

#define MAXSEQLEN 65536

#define FALSE 0
#define TRUE 1

#define SQR(x) ((x)*(x))
#define MIN(x,y) (((x)<(y))?(x):(y))
#define MAX(x,y) (((x)>(y))?(x):(y))

const char *rescodes = "ARNDCQEGHILKMFPSTWYVBZX";

/*  BLOSUM 62 */
const short           aamat[23][23] =
{
    {4, -1, -2, -2, 0, -1, -1, 0, -2, -1, -1, -1, -1, -2, -1, 1, 0, -3, -2, 0, -2, -1, 0},
    {-1, 5, 0, -2, -3, 1, 0, -2, 0, -3, -2, 2, -1, -3, -2, -1, -1, -3, -2, -3, -1, 0, -1},
    {-2, 0, 6, 1, -3, 0, 0, 0, 1, -3, -3, 0, -2, -3, -2, 1, 0, -4, -2, -3, 3, 0, -1},
    {-2, -2, 1, 6, -3, 0, 2, -1, -1, -3, -4, -1, -3, -3, -1, 0, -1, -4,
     -3, -3, 4, 1, -1},
    {0, -3, -3, -3,10, -3, -4, -3, -3, -1, -1, -3, -1, -2, -3, -1, -1, -2,
     -2, -1, -3, -3, -2},
    {-1, 1, 0, 0, -3, 5, 2, -2, 0, -3, -2, 1, 0, -3, -1, 0, -1, -2,
     -1, -2, 0, 3, -1},
    {-1, 0, 0, 2, -4, 2, 5, -2, 0, -3, -3, 1, -2, -3, -1, 0, -1, -3,
     -2, -2, 1, 4, -1},
    {0, -2, 0, -1, -3, -2, -2, 6, -2, -4, -4, -2, -3, -3, -2, 0, -2, -2,
     -3, -3, -1, -2, -1},
    {-2, 0, 1, -1, -3, 0, 0, -2, 8, -3, -3, -1, -2, -1, -2, -1, -2, -2,
     2, -3, 0, 0, -1},
    {-1, -3, -3, -3, -1, -3, -3, -4, -3, 4, 2, -3, 1, 0, -3, -2, -1, -3,
     -1, 3, -3, -3, -1},
    {-1, -2, -3, -4, -1, -2, -3, -4, -3, 2, 4, -2, 2, 0, -3, -2, -1, -2,
     -1, 1, -4, -3, -1},
    {-1, 2, 0, -1, -3, 1, 1, -2, -1, -3, -2, 5, -1, -3, -1, 0, -1, -3,
     -2, -2, 0, 1, -1},
    {-1, -1, -2, -3, -1, 0, -2, -3, -2, 1, 2, -1, 5, 0, -2, -1, -1, -1,
     -1, 1, -3, -1, -1},
    {-2, -3, -3, -3, -2, -3, -3, -3, -1, 0, 0, -3, 0, 6, -4, -2, -2, 1,
     3, -1, -3, -3, -1},
    {-1, -2, -2, -1, -3, -1, -1, -2, -2, -3, -3, -1, -2, -4, 7, -1, -1, -4,
     -3, -2, -2, -1, -2},
    {1, -1, 1, 0, -1, 0, 0, 0, -1, -2, -2, 0, -1, -2, -1, 4, 1, -3,
     -2, -2, 0, 0, 0},
    {0, -1, 0, -1, -1, -1, -1, -2, -2, -1, -1, -1, -1, -2, -1, 1, 5, -2,
     -2, 0, -1, -1, 0},
    {-3, -3, -4, -4, -2, -2, -3, -2, -2, -3, -2, -3, -1, 1, -4, -3, -2, 11,
     2, -3, -4, -3, -2},
    {-2, -2, -2, -3, -2, -1, -2, -3, 2, -1, -1, -2, -1, 3, -3, -2, -2, 2,
     7, -1, -3, -2, -1},
    {0, -3, -3, -3, -1, -2, -2, -3, -3, 3, 1, -2, 1, -1, -2, -2, 0, -3,
     -1, 4, -3, -2, -1},
    {-2, -1, 3, 4, -3, 0, 1, -1, 0, -3, -4, 0, -3, -3, -2, 0, -1, -4,
     -3, -3, 4, 1, -1},
    {-1, 0, 0, 1, -3, 3, 4, -2, 0, -3, -3, 1, -1, -3, -1, 0, -1, -3,
     -2, -2, 1, 4, -1},
    {0, -1, -1, -1, -2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -2, 0, 0, -2,
     -1, -1, -1, -1, 4}
};

/* Dump a rude message to standard error and exit */
void
  fail(char *errstr)
{
    fprintf(stderr, "\n*** %s\n\n", errstr);
    exit(-1);
}

/* Convert AA letter to numeric code (0-22) */
int
  aanum(int ch)
{
    static int aacvs[] =
    {
	999, 0, 20, 4, 3, 6, 13, 7, 8, 9, 22, 11, 10, 12, 2,
	22, 14, 5, 1, 15, 16, 22, 19, 17, 22, 18, 21
    };

    return (isalpha(ch) ? aacvs[ch & 31] : 22);
}

/* This routine will read in one sequence from a database file. The
   sequence can be in any of the supported formats. Returns length
   of sequence.
*/
int
  getseq(char *dbname, char *dseq, FILE * lfil)
{
    int i, j, len;
    short badln, fformat;
    enum
    {
	UNKNOWN, EMBL, FASTA, OWL, GCG
    };
    char buf[MAXSEQLEN], split;
    int offset;

    offset = j = 0;

    if (!fgets(buf, MAXSEQLEN, lfil))
	return (-1);
    if (strstr(buf, "of:") != NULL && strstr(buf, "check:") != NULL)
	fformat = GCG;
    else if (strncmp(buf, "ID   ", 5) == 0)
	fformat = EMBL;
    else if (buf[0] == '>' && (buf[1] == '>' || buf[3] == ';'))
	fformat = OWL;
    else if (buf[0] == '>')
	fformat = FASTA;
    else
    {
	fprintf(stderr, "WARNING: Attempting to interpret input file with unknown format");
	fformat = UNKNOWN;
    }

    switch (fformat)
    { 
    case GCG:
	sscanf(strstr(buf, "of:")+3, "%s", dbname);
	while (strstr(buf, "..") == NULL)
	    fgets(buf, MAXSEQLEN, lfil);
	fgets(buf, MAXSEQLEN, lfil);
	break;
	
    case EMBL:
	strncpy(dbname, buf + 5, 70);
	while (buf[0] != ' ')
	    fgets(buf, MAXSEQLEN, lfil);
	break;
	
    case OWL:
	fgets(buf, MAXSEQLEN, lfil);
	strncpy(dbname, buf, 70);
	fgets(buf, MAXSEQLEN, lfil);
	break;
	
    case FASTA:
	strncpy(dbname, buf + 1, 70);
	fgets(buf, MAXSEQLEN, lfil);
	break;
	
    default:
	/* Try to find a line which looks like a protein sequence */
	do
	{
	    badln = (strpbrk(buf, "JjOoUu<>#$%&@") != NULL);
	    if (badln && !fgets(buf, MAXSEQLEN, lfil))
		return (-1);
	}
	while (badln);
	strcpy(dbname, "<NO NAME>");
	break;
    }

    if (dbname[(len = strlen(dbname)) - 1] == '\n')
	dbname[--len] = '\0';
    if (len >= 70)
	dbname[70] = '\0';

    for (;;)
    {
	if (!strncmp(buf, "//", 2))
	    break;
	len = strlen(buf);
	for (i = offset; i < len && j < MAXSEQLEN; i++)
	{
	    split = islower(buf[i]) ? toupper(buf[i]) : buf[i];
	    if (split == '@' || (fformat == OWL && split == '*'))
	    {
		dseq[j] = '\0';
		while (fgets(buf, MAXSEQLEN, lfil));
		return (j);
	    }
	    if (isalpha(split))
		dseq[j++] = split;
	    else if (buf[i] == '\n')
		break;
	}
	if (!fgets(buf, MAXSEQLEN, lfil))
	    break;
    }

    if (j == MAXSEQLEN)
	printf("\nWARNING: sequence %s over %d long; truncated!\n",
	       dbname, MAXSEQLEN);

    dseq[j] = '\0';
    return (j);
}

int main(int argc, char **argv)
{
    int i, j, seqlen=0;
    char desc[65536], seq[MAXSEQLEN], buf[65536], *p;
    char *ncbicodes = "XAXCDEFGHIKLMNPQRSTVWXYXXX";
    FILE *ifp;

    if (argc != 2)
	fail("Usage: seq2psi seq-file");

    ifp = fopen(argv[1], "r");
    if (!ifp)
	fail("Unable to open sequence file!");

    seqlen = getseq(desc, seq, ifp);

    if (seqlen < 5 || seqlen >= MAXSEQLEN)
	fail("Sequence length error!");

    printf("%d\n", seqlen);

    for (i=0; i<seqlen; i++)
      putchar(seq[i]);

    printf("\n0\n0\n0\n0\n0\n0\n0\n0\n0\n0\n0\n0\n");

    for (i=0; i<seqlen; i++)
    {
	for (j=0; j<26; j++)
	  if (ncbicodes[j] != 'X')
	    printf("%d  ", aamat[aanum(seq[i])][aanum(ncbicodes[j])]*100);
	  else
	    printf("-32768  ");
	putchar('\n');
    }

    return 0;
}