/* program:  すべてのアミノ酸配列同士の seq. identity 計算
             注）seq. id. は切捨て
                 長さが 1500 以上のものは（入力するが）無視
   file:     seqid.c
   date:     2003/7/3 by O.W
*/

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

#define MAXINLENG 2500  /* 入力配列長最大値 */
#define MAXLENG   1500  /* 配列長最大値 */
#define MAXDATA   1600  /* 入力総タンパク質数上限 */
#define ASIZE 23        /* アミノ酸の種類（種類は 22 個，その他の分類が 1 個）*/
#define INFTY -70000    /* 無限値（ありえないスコア）*/
#define GOPEN -10       /* ギャップを開くときのペナルティ */
#define GEXT  -1        /* ギャップを続けるときのペナルティ */

/* スコア表 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 score(int, int);

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 sid, i;

  FILE *fp;
  char instr[MAXINLENG];

  /* データ列の読み込み */
  fp = fopen("proteins.txt", "r");
  if(!fp) {
    printf("File open error!\n");
    exit(1);
  }
  num = 0;  tleng = 0;
  while(1) {
    fgets(instr, MAXINLENG + 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'];
    }
    if(i > MAXLENG)
      LENG[num] = -i;     // 長さを負にすることで以後無視する
    else
      LENG[num] = i;
    tleng = tleng + i;
  }
  fclose(fp);
  mleng = tleng / num;

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

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

  return 0;
}

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

  int i, j;

  /* 探索打ち切り */
  /* 解説：一方の列が終わったらペナルティ 0 で打ち切ってよい */
  for(i = 1; i <= m; i++)
    STB[i][n][0] = TBL[ALL[k1][i]][ALL[k2][n]];
  for(j = 1; j <= n; j++)
    STB[m][j][0] = TBL[ALL[k1][m]][ALL[k2][j]];
  /* 探索打ち切り */
  /* 解説：ギャップ途中で一方の列が終わることは許されない */
  for(i = 1; i <= m; i++)
    STB[i][n][1] = TBL[ALL[k1][i]][ALL[k2][n]];
  for(j = 1; j <= n; j++)
    STB[m][j][2] = TBL[ALL[k1][m]][ALL[k2][j]];

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

      STB[i][j][1] = max(STB[i+1][j+1][0] + TBL[ch][ALL[k2][j]],
			 STB[ i ][j+1][1] + GEXT);

      STB[i][j][2] = max(STB[i+1][j+1][0] + TBL[ch][ALL[k2][j]],
			 STB[i+1][j][2] + GEXT);
    }
  }

  return;
}
      
int score(int m, int n) {
  int i, j, s1, s2;

  /* 最適スタート地点を探す */
  s1 = INFTY;
  for(j = 1; j <= n; j++) s1 = max(s1, STB[1][j][0]);

  s2 = INFTY;
  for(i = 1; i <= m; i++) s2 = max(s2, STB[i][1][0]);

  return max(s1, s2);
}