/* program:  すべてのアミノ酸配列同士の seq. identity 計算
             注）seq. id. は切捨て
   file:     seqid.c
   date:     2003/7/23 by O.W
*/

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

#define MAXLENG 3000    /* 入力配列長最大値 */
#define MAXDATA 1600    /* 入力総タンパク質数上限 */
#define ASIZE 23        /* アミノ酸の種類（種類は 22 個，他の分類が 1 個）*/
#define INFTY -70000000 /* 無限値（ありえないスコア）NEW */
#define SHIFT  10000    /* スコアの値をシフトするための係数 NEW */
#define GOPEN -100000   /* ギャップを開くときのペナルティ * SHIFT NEW */
#define GEXT  -10000    /* ギャップを続けるときのペナルティ * SHIFT NEW */

/* スコア表 BLOSUM-62 */
int TBL[ASIZE+1][ASIZE+1] = {
  { -9,-9,-9,-9,-9,-9,-9,-9,-9,-9,-9,-9,-9,-9,-9,-9,-9,-9,-9,-9,-9,-9,-9,-9 },
  { -9, 4,-1,-2,-2, 0,-1,-1, 0,-2,-1,-1,-1,-1,-2,-1, 1, 0,-3,-2, 0,-2,-1, 0 },
  { -9,-1, 5, 0,-2,-3, 1, 0,-2, 0,-3,-2, 2,-1,-3,-2,-1,-1,-3,-2,-3,-1, 0,-1 },
  { -9,-2, 0, 6, 1,-3, 0, 0, 0, 1,-3,-3, 0,-2,-3,-2, 1, 0,-4,-2,-3, 3, 0,-1 },
  { -9,-2,-2, 1, 6,-3, 0, 2,-1,-1,-3,-4,-1,-3,-3,-1, 0,-1,-4,-3,-3, 4, 1,-1 },
  { -9, 0,-3,-3,-3, 9,-3,-4,-3,-3,-1,-1,-3,-1,-2,-3,-1,-1,-2,-2,-1,-3,-3,-2 },
  { -9,-1, 1, 0, 0,-3, 5, 2,-2, 0,-3,-2, 1, 0,-3,-1, 0,-1,-2,-1,-2, 0, 3,-1 },
  { -9,-1, 0, 0, 2,-4, 2, 5,-2, 0,-3,-3, 1,-2,-3,-1, 0,-1,-3,-2,-2, 1, 4,-1 },
  { -9, 0,-2, 0,-1,-3,-2,-2, 6,-2,-4,-4,-2,-3,-3,-2, 0,-2,-2,-3,-3,-1,-2,-1 },
  { -9,-2, 0, 1,-1,-3, 0, 0,-2, 8,-3,-3,-1,-2,-1,-2,-1,-2,-2, 2,-3, 0, 0,-1 },
  { -9,-1,-3,-3,-3,-1,-3,-3,-4,-3, 4, 2,-3, 1, 0,-3,-2,-1,-3,-1, 3,-3,-3,-1 },
  { -9,-1,-2,-3,-4,-1,-2,-3,-4,-3, 2, 4,-2, 2, 0,-3,-2,-1,-2,-1, 1,-4,-3,-1 },
  { -9,-1, 2, 0,-1,-3, 1, 1,-2,-1,-3,-2, 5,-1,-3,-1, 0,-1,-3,-2,-2, 0, 1,-1 },
  { -9,-1,-1,-2,-3,-1, 0,-2,-3,-2, 1, 2,-1, 5, 0,-2,-1,-1,-1,-1, 1,-3,-1,-1 },
  { -9,-2,-3,-3,-3,-2,-3,-3,-3,-1, 0, 0,-3, 0, 6,-4,-2,-2, 1, 3,-1,-3,-3,-1 },
  { -9,-1,-2,-2,-1,-3,-1,-1,-2,-2,-3,-3,-1,-2,-4, 7,-1,-1,-4,-3,-2,-2,-1,-2 },
  { -9, 1,-1, 1, 0,-1, 0, 0, 0,-1,-2,-2, 0,-1,-2,-1, 4, 1,-3,-2,-2, 0, 0, 0 },
  { -9, 0,-1, 0,-1,-1,-1,-1,-2,-2,-1,-1,-1,-1,-2,-1, 1, 5,-2,-2, 0,-1,-1, 0 },
  { -9,-3,-3,-4,-4,-2,-2,-3,-2,-2,-3,-2,-3,-1, 1,-4,-3,-2,11, 2,-3,-4,-3,-2 },
  { -9,-2,-2,-2,-3,-2,-1,-2,-3, 2,-1,-1,-2,-1, 3,-3,-2,-2, 2, 7,-1,-3,-2,-1 },
  { -9, 0,-3,-3,-3,-1,-2,-2,-3,-3, 3, 1,-2, 1,-1,-2,-2, 0,-3,-1, 4,-3,-2,-1 },
  { -9,-2,-1, 3, 4,-3, 0, 1,-1, 0,-3,-4, 0,-3,-3,-2, 0,-1,-4,-3,-3, 4, 1,-1 },
  { -9,-1, 0, 0, 1,-3, 3, 4,-2, 0,-3,-3, 1,-1,-3,-1, 0,-1,-3,-2,-2, 1, 4,-1 },
  { -9, 0,-1,-1,-1,-2,-1,-1,-1,-1,-1,-1,-1,-1,-1,-2, 0, 0,-2,-1,-1,-1,-1,-1 }
};
/* 補足：
  1. BLAST http://www.ncbi.nlm.nih.gov/blast/html/sub_matrix.html より
  2. アルファベット（アミノ酸記号）とインデックスの対応は以下の表の通り
*/
/* アルファベット表（X は「その他」用）*/
int C2N[26] = { /* char --> index */
   1,21, 5, 4, 7,14, 8, 9,10, 0,12,11,13, 3, 0,15, 6, 2,16,17, 0,20,18,23,19,22
/* A, B, C, D, E, F, G, H, I, J, K, L, M, N, O, P, Q, R, S, T, U, V, W, X, Y, Z */
};
int N2C[24] = { /* index --> char */
   '*', 'A', 'R', 'N', 'D', 'C', 'Q', 'E', 'G', 'H', 'I', 'L',
   'K', 'M', 'F', 'P', 'S', 'T', 'W', 'Y', 'V', 'B', 'Z', 'X'
};

/* 大域変数 */
int ALL[MAXDATA+1][MAXLENG+1];     /* ALL[k] = k 本目のアミノ酸配列 */
int LENG[MAXDATA+1];               /* LENG[k] = ALL[k] の長さ */
int STB[MAXLENG+1][MAXLENG+1][3];  /* スコア計算のための表 */

void maketb(int, int, int, int);
int idnum(int, int, int, int);

int match(int c1, int c2) { return (c1 != 23) & (c1 == c2); }
int max(x, y)     { return x > y ? x : y; }
int max3(x, y, z) { return max(max(x, y), z); }
int min(x, y)     { return x < y ? x : y; }

int main(int argc, char *argv[]) {

  int leng, tleng, mleng;
  int num;
  int m, n, k1, k2;
  int i, j, sid;

  FILE *fp;
  char instr[MAXLENG];

  /* データ列の読み込み */
  fp = fopen("proteins-smallexample.txt", "r");
  //fp = fopen("proteinsA.txt", "r");
  //fp = fopen("proteins.txt", "r");
  if(!fp) {
    printf("File open error!\n");
    exit(1);
  }
  num = 0;  tleng = 0;
  while(1) {
    fgets(instr, MAXLENG + 1, fp);
    if(instr[0] == '@') break;
    num++;
    leng = strlen(instr);
    for(i = 0; i < leng; i++) {
      if(instr[i] == '%') break;
      ALL[num][i+1] = C2N[instr[i] - 'A'];
    }
    LENG[num] = i;
    tleng = tleng + i;
  }
  fclose(fp);
  mleng = tleng / num;

  printf("number of seq.s = %d,  average length = %d\n", num, mleng);

  /* TBL の変換（SHIFT 倍する）NEW */
  for(i = 0; i <= ASIZE; i++)
    for(j = 0; j <= ASIZE; j++)
      TBL[i][j] = TBL[i][j] * SHIFT;

  /* alignment, seq. id.  の計算 */
  for(k1 = 1; k1 <= num; k1++) {
    m = LENG[k1];
    printf("%3d, %4d:", k1, m);
    for(k2 = k1 + 1; k2 <= num; k2++) {
      n = LENG[k2];
      maketb(k1, k2, m, n);
      sid = idnum(k1, k2, m, n) * 100 / min(m, n);
      printf(" %2d,", sid);
    }
    printf("\n");
  }

  return 0;
}

void maketb(int k1, int k2, int m, int n) {

  int i, j, ch, ch1, ch2, t;

  /* 探索打ち切りのための境界条件 NEW */
  /* 解説：
     1. ギャップ途中で一方の列が終わることは許されない ==> 最後は必ず比較
     2. 通常の時に一方の列が終わったらペナルティ 0 で打ち切ってよい
  */
  ch = ALL[k2][n];
  for(i = 1; i <= m; i++)
    STB[i][n][0] = STB[i][n][1] = TBL[ALL[k1][i]][ch] + match(ALL[k1][i],ch);
  ch = ALL[k1][m];
  for(j = 1; j <= n; j++)
    STB[m][j][0] = STB[m][j][2] = TBL[ch][ALL[k2][j]] + match(ALL[k2][j],ch);

  /* STB の計算 */
  for(i = m - 1; i >= 1; i--) {
    ch1 = ALL[k1][i];
    for(j = n - 1; j >= 1; j--) {
      ch2 = ALL[k2][j];
      t = TBL[ch1][ch2];
      t = t + match(ch1, ch2);  // NEW
      STB[i][j][0] = max3(STB[i+1][j+1][0] + t,
			  STB[ i ][j+1][1] + GOPEN,
			  STB[i+1][ j ][2] + GOPEN);

      STB[i][j][1] = max(STB[i+1][j+1][0] + t,
			 STB[ i ][j+1][1] + GEXT);

      STB[i][j][2] = max(STB[i+1][j+1][0] + t,
			 STB[i+1][ j ][2] + GEXT);
    }
  }

  return;
}

int idnum(int k1, int k2, int m, int n) {

  int i, j;
  int maxs, ans;

  /* 最大 mixed score を求める NEW */
  maxs = INFTY;
  for(j = 1; j <= n; j++) maxs = max(maxs, STB[1][j][0]);
  for(i = 1; i <= m; i++) maxs = max(maxs, STB[i][1][0]);

  /* 一致個数 idnum を取り出す NEW */
  if(maxs < 0) ans = maxs + (int)(-(float)maxs/(float)SHIFT + 0.5) * SHIFT;
  else         ans = maxs % SHIFT;
  
  return ans;
}