/*****************************************************************
 * 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.
 *****************************************************************/

/* clustal.c
 * SRE, Sun Jun  6 17:50:45 1999 [bus from Madison, 1999 worm mtg]
 * 
 * Import/export of ClustalV/W multiple sequence alignment
 * formatted files. Derivative of msf.c; MSF is a pretty
 * generic interleaved format.   
 * 
 * RCS $Id: clustal.c 228 2011-03-29 14:05:27Z dave $ (Original squid RCS Id: clustal.c,v 1.1 1999/07/15 22:26:53 eddy Exp)
 */

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <ctype.h>
#include "squid.h"
#include "msa.h"

#ifdef CLUSTALO
/* needed for PACKAGE_VERSION */
#include "../config.h"

#ifndef min
	#define min( a, b ) ( ((a) < (b)) ? (a) : (b) )
#endif

/*These are all the positively scoring groups that occur in the Gonnet Pam250
matrix. There are strong and weak groups, defined as strong score >0.5 and
weak score =<0.5. Strong matching columns to be assigned ':' and weak matches
assigned '.' in the clustal output format.
amino_strong = res_cat1
amino_weak = res_cat2
*/

char *amino_strong[] = {"STA", "NEQK", "NHQK", "NDEQ", "QHRK", "MILV",
    "MILF", "HY", "FYW", NULL};

char *amino_weak[] = {"CSA", "ATV", "SAG", "STNK", "STPA", "SGND",
    "SNDEQK", "NDEQHK", "NEQHRK", "FVLIM", "HFY", NULL};

#endif

#ifdef TESTDRIVE_CLUSTAL
/*****************************************************************
 * msf.c test driver: 
 * cc -DTESTDRIVE_CLUSTAL -g -O2 -Wall -o test clustal.c msa.c gki.c sqerror.c sre_string.c file.c hsregex.c sre_math.c sre_ctype.c -lm
 * 
 */
int
main(int argc, char **argv)
{
  MSAFILE *afp;
  MSA     *msa;
  char    *file;
  
  file = argv[1];

  if ((afp = MSAFileOpen(file, MSAFILE_CLUSTAL, NULL)) == NULL)
    Die("Couldn't open %s\n", file);

  while ((msa = ReadClustal(afp)) != NULL)
    {
      WriteClustal(stdout, msa);
      MSAFree(msa); 
    }
  
  MSAFileClose(afp);
  exit(0);
}
/******************************************************************/
#endif /* testdrive_clustal */


/* Function: ReadClustal()
 * Date:     SRE, Sun Jun  6 17:53:49 1999 [bus from Madison, 1999 worm mtg]
 *
 * Purpose:  Parse an alignment read from an open Clustal format
 *           alignment file. Clustal is a single-alignment format.
 *           Return the alignment, or NULL if we have no data.
 *           
 * Args:     afp  - open alignment file
 *
 * Returns:  MSA * - an alignment object
 *                   caller responsible for an MSAFree()
 *           NULL if no more alignments
 *
 * Diagnostics: 
 *           Will Die() here with a (potentially) useful message
 *           if a parsing error occurs.
 */
MSA *
ReadClustal(MSAFILE *afp)
{
  MSA    *msa;
  char   *s;
  int     slen;
  int     sqidx;
  char   *name;
  char   *seq;
  char   *s2;

  if (feof(afp->f)) return NULL;

  /* Skip until we see the CLUSTAL header
   */
  while ((s = MSAFileGetLine(afp)) != NULL)
    {
      if (strncmp(s, "CLUSTAL", 7) == 0 &&
	  strstr(s, "multiple sequence alignment") != NULL)
	break;
    }
  if (s == NULL) return NULL;

  msa = MSAAlloc(10, 0);

  /* Now we're in the sequence section. 
   * As discussed above, if we haven't seen a sequence name, then we
   * don't include the sequence in the alignment.
   * Watch out for conservation markup lines that contain *.: chars
   */
  while ((s = MSAFileGetLine(afp)) != NULL) 
    {
      if ((name = sre_strtok(&s, WHITESPACE, NULL))  == NULL) continue;
      if ((seq  = sre_strtok(&s, WHITESPACE, &slen)) == NULL) continue;
      s2 = sre_strtok(&s, "\n", NULL);

      /* The test for a conservation markup line
       */
      if (strpbrk(name, ".*:") != NULL && strpbrk(seq, ".*:") != NULL)
	continue;
      if (s2 != NULL)
	Die("Parse failed at line %d, file %s: possibly using spaces as gaps",
	    afp->linenumber, afp->fname);
  
      /* It's not blank, and it's not a coord line: must be sequence
       */
      sqidx = MSAGetSeqidx(msa, name, msa->lastidx+1);
      msa->lastidx = sqidx;
      msa->sqlen[sqidx] = sre_strcat(&(msa->aseq[sqidx]), msa->sqlen[sqidx], seq, slen); 
    }

  MSAVerifyParse(msa);		/* verifies, and also sets alen and wgt. */
  return msa;
}


/* Function: WriteClustal()
 * Date:     SRE, Sun Jun  6 18:12:47 1999 [bus from Madison, worm mtg 1999]
 *
 * Purpose:  Write an alignment in Clustal format to an open file.
 *
 * Args:     fp    - file that's open for writing.
 *           msa   - alignment to write. 
 *
 * Returns:  (void)
 */
void
WriteClustal(FILE *fp, MSA *msa)
{
  int    idx;			/* counter for sequences         */
  int    len;			/* tmp variable for name lengths */
  int    namelen;		/* maximum name length used      */
  int    pos;			/* position counter              */
  char   buf[80];	        /* buffer for writing seq        */
  int    cpl = 60;		/* char per line (< 64)          */

  /* consensus line stuff */
  int subpos;
  char first;
  int bail;
  int strong_bins[9];
  int weak_bins[11];
  int cons;
  int bin;

  /* calculate max namelen used */
  namelen = 0;
  for (idx = 0; idx < msa->nseq; idx++)
    if ((len = strlen(msa->sqname[idx])) > namelen) 
      namelen = len;

#ifdef CLUSTALO
  fprintf(fp, "CLUSTAL O(%s) multiple sequence alignment\n", PACKAGE_VERSION);
#else
  fprintf(fp, "CLUSTAL W(1.5) multiple sequence alignment\n");
#endif
  
  /*****************************************************
   * Write the sequences
   *****************************************************/

#ifdef CLUSTALO
    fprintf(fp, "\n");	/* original had two blank lines */
#endif

  for (pos = 0; pos < msa->alen; pos += cpl)
    {
      fprintf(fp, "\n");	/* Blank line between sequence blocks */
      for (idx = 0; idx < msa->nseq; idx++)
      {
        strncpy(buf, msa->aseq[idx] + pos, cpl);
	    buf[cpl] = '\0';
#ifdef CLUSTALO
	    fprintf(fp, "%-*s %s\n", namelen+5, msa->sqname[idx], buf);
#else
	    fprintf(fp, "%*s %s\n", namelen, msa->sqname[idx], buf);
#endif
	  }
#ifdef CLUSTALO
      /* do consensus dots */

      /* print namelen+5 spaces */
      for(subpos = 0; subpos <= namelen+5; subpos++)
        fprintf(fp, " ");

      for(subpos = pos; subpos < min(pos + cpl, msa->alen); subpos++)
      {
          /* see if 100% conservation */
          first = msa->aseq[0][subpos];
          bail = 0;
          for (idx = 1; idx < msa->nseq; idx++)
          {
            if(msa->aseq[idx][subpos] != first)
            {
              bail = 1;
              break;
            }
          }
          if(!bail)
            fprintf(fp, "*");
          else
          {
            /* if not then check strong */
            for(bin = 0; bin < 9; bin++)
              strong_bins[bin] = 0; /* clear the bins */

            for(idx = 0; idx < msa->nseq; idx++)
            {
              switch(msa->aseq[idx][subpos])
              {
                case 'S': strong_bins[0]++; break;
                case 'T': strong_bins[0]++; break;
                case 'A': strong_bins[0]++; break;
                case 'N': strong_bins[1]++; strong_bins[2]++; strong_bins[3]++; break;
                case 'E': strong_bins[1]++; strong_bins[3]++; break;
                case 'Q': strong_bins[1]++; strong_bins[2]++; strong_bins[3]++; strong_bins[4]++; break;
                case 'K': strong_bins[1]++; strong_bins[2]++; strong_bins[4]++; break;
                case 'D': strong_bins[3]++; break;
                case 'R': strong_bins[4]++; break;
                case 'H': strong_bins[4]++; strong_bins[7]++; break;
                case 'M': strong_bins[5]++; strong_bins[6]++; break;
                case 'I': strong_bins[5]++; strong_bins[6]++; break;
                case 'L': strong_bins[5]++; strong_bins[6]++; break;
                case 'V': strong_bins[5]++; break;
                case 'F': strong_bins[6]++; strong_bins[8]++; break;
                case 'Y': strong_bins[7]++; strong_bins[8]++; break;
                case 'W': strong_bins[8]++; break;
              }
            }
            bail = 0;
            for(bin = 0; bin < 9; bin++)
              if(strong_bins[bin] == msa->nseq)
              {
                  bail = 1;
                  break;
              }
            if(bail)
              fprintf(fp, ":");
            else
            {
              /* check weak */
              for(bin = 0; bin < 11; bin++)
                weak_bins[bin] = 0; /* clear the bins */

              for(idx = 0; idx < msa->nseq; idx++)
              {
                switch(msa->aseq[idx][subpos])
                {
                  case 'C': weak_bins[0]++; break;
                  case 'S': weak_bins[0]++; weak_bins[2]++; weak_bins[3]++; weak_bins[4]++; weak_bins[5]++; weak_bins[6]++; break;
                  case 'A': weak_bins[0]++; weak_bins[1]++; weak_bins[2]++; weak_bins[4]++; break;
                  case 'T': weak_bins[1]++; weak_bins[3]++; weak_bins[4]++; break;
                  case 'V': weak_bins[1]++; weak_bins[9]++; break;
                  case 'G': weak_bins[2]++; break;
                  case 'N': weak_bins[3]++; weak_bins[5]++; weak_bins[6]++; weak_bins[7]++; weak_bins[8]++; break;
                  case 'K': weak_bins[3]++; weak_bins[6]++; weak_bins[7]++; weak_bins[8]++; break;
                  case 'D': weak_bins[5]++; weak_bins[6]++; weak_bins[7]++; break;
                  case 'E': weak_bins[6]++; weak_bins[7]++; weak_bins[8]++; break;
                  case 'Q': weak_bins[6]++; weak_bins[7]++; weak_bins[8]++; break;
                  case 'H': weak_bins[7]++; weak_bins[8]++; weak_bins[10]++; break;
                  case 'R': weak_bins[8]++; break;
                  case 'F': weak_bins[9]++; weak_bins[10]++; break;
                  case 'L': weak_bins[9]++; break;
                  case 'I': weak_bins[9]++; break;
                  case 'M': weak_bins[9]++; break;
                  case 'Y': weak_bins[10]++; break;
                }
              }
              bail = 0;
              for(bin = 0; bin < 11; bin++)
                if(weak_bins[bin] == msa->nseq)
                {
                    bail = 1;
                    break;
                }
              if(bail)
                fprintf(fp, ".");
              else
                fprintf(fp, " ");
            }
          }
      }
      fprintf(fp,"\n");
#endif
    }

  return;
}