binaries/src/GLProbs-1.0/MSAReadMatrix.cpp

   1 /***********************************************
   2  * # Copyright 2009-2010. Liu Yongchao
   3  * # Contact: Liu Yongchao, School of Computer Engineering,
   4  * #                     Nanyang Technological University.
   5  * # Emails:     liuy0039@ntu.edu.sg; nkcslyc@hotmail.com
   6  * #
   7  * # GPL version 3.0 applies.
   8  * #
   9  * ************************************************/
  10
  11 #include <string.h>
  12 #include <stdio.h>
  13 #include <stdlib.h>
  14 #include <math.h>
  15 #include "MSAReadMatrix.h"
  16
  17 #define TRACE 0
  18
  19 ////////////////////////////////////////////////////////////
  20 // extern variables for scoring matrix data
  21 ////////////////////////////////////////////////////////////
  22 extern float g_gap_open1, g_gap_open2, g_gap_ext1, g_gap_ext2;
  23 extern char *aminos, *bases, matrixtype[20];
  24 extern int subst_index[26];
  25
  26 extern double sub_matrix[26][26];
  27 extern double normalized_matrix[26][26];
  28
  29 extern float TEMPERATURE;
  30 extern int MATRIXTYPE;
  31
  32 extern float GAPOPEN;
  33 extern float GAPEXT;
  34
  35 typedef struct {
  36         char input[30];
  37         int matrix;
  38         int N;
  39         float T;
  40         float beta;
  41         char opt;                       //can be 'P' or 'M'
  42         float gapopen;
  43         float gapext;
  44 } argument_decl;
  45
  46 //argument support
  47 extern argument_decl argument;
  48
  49 /////////////////////////////////////////////////////////
  50 //sets substitution matrix type
  51 ////////////////////////////////////////////////////////
  52 void setmatrixtype(int le) {
  53         switch (le) {
  54         case 160:
  55                 strcpy(matrixtype, "gonnet_160");
  56                 break;
  57         case 4:
  58                 strcpy(matrixtype, "nuc_simple");
  59                 break;
  60         default:
  61                 strcpy(matrixtype, "CUSTOM");
  62                 break;
  63
  64         };
  65
  66 }
  67
  68 ///////////////////////////////////////////////////////////////////
  69 //sets matrix flag
  70 ///////////////////////////////////////////////////////////////////
  71 inline int matrixtype_to_int() {
  72
  73         if (!strcmp(matrixtype, "nuc_simple"))
  74                 return 4;
  75         else if (!strcmp(matrixtype, "gonnet_160"))
  76                 return 160;
  77         else
  78                 return 1000;
  79
  80 }
  81
  82 /////////////////////////////////////////////////////////////////
  83 //
  84 // Can read any scoring matrix as long as it is defined in Matrix.h
  85 // AND it is a lower triangular
  86 // AND the order of amino acids/bases is mentioned
  87 /////////////////////////////////////////////////////////////////
  88
  89 inline void read_matrix(score_matrix matrx) {
  90         int i, j, basecount, position = 0;
  91
  92         bases = (char *) matrx.monomers;
  93
  94         basecount = strlen(bases);
  95
  96         for (i = 0; i < basecount; i++)
  97                 subst_index[i] = -1;
  98
  99         for (i = 0; i < basecount; i++)
 100                 subst_index[bases[i] - 'A'] = i;
 101
 102         if (TRACE == 1)
 103                 printf("\nbases read: %d\n", basecount);
 104
 105         for (i = 0; i < basecount; i++)
 106                 for (j = 0; j <= i; j++) {
 107
 108                         double value = exp(argument.beta * matrx.matrix[position++]);
 109                         sub_matrix[i][j] = value;
 110                         sub_matrix[j][i] = value;
 111                 }
 112
 113         if (TRACE)
 114                 for (i = 0; i < basecount; i++) {
 115                         for (j = 0; j < basecount; j++)
 116                                 printf(" %g ", sub_matrix[i][j]);
 117                         printf("\n");
 118                 }
 119
 120 }
 121
 122 /////////////////////////////////////////////////////////////////
 123 // read normalized residue exchange matrix
 124 // compute sequence similarity
 125 // add by YE Yongtao
 126 /////////////////////////////////////////////////////////////////
 127
 128 inline void read_normalized_matrix(score_matrix matrx) {
 129         int i, j, basecount, position = 0;
 130
 131         bases = (char *) matrx.monomers;
 132
 133         basecount = strlen(bases);
 134
 135         for (i = 0; i < basecount; i++)
 136                 subst_index[i] = -1;
 137
 138         for (i = 0; i < basecount; i++)
 139                 subst_index[bases[i] - 'A'] = i;
 140
 141         if (TRACE == 1)
 142                 printf("\nbases read: %d\n", basecount);
 143
 144         for (i = 0; i < basecount; i++)
 145                 for (j = 0; j <= i; j++) {
 146
 147                         double value =  matrx.matrix[position++];
 148                         normalized_matrix[i][j] = value;
 149                         normalized_matrix[j][i] = value;
 150                 }
 151
 152         if (TRACE)
 153                 for (i = 0; i < basecount; i++) {
 154                         for (j = 0; j < basecount; j++)
 155                                 printf(" %g ", normalized_matrix[i][j]);
 156                         printf("\n");
 157                 }
 158
 159 }
 160 //////////////////////////////////////////////////////////////////////////////////
 161 //intialize the arguments (default values)
 162 //////////////////////////////////////////////////////////////////////////////////
 163 void init_arguments() {
 164         float gap_open = 0, gap_ext = 0;
 165         int le;
 166
 167         le = matrixtype_to_int();
 168
 169         argument.N = 1;
 170         strcpy(argument.input, "tempin");
 171         argument.matrix = le;
 172         argument.gapopen = GAPOPEN;
 173         argument.gapext = GAPEXT;
 174         argument.T = TEMPERATURE;
 175         argument.beta = 1.0 / TEMPERATURE;
 176         argument.opt = 'P';
 177
 178         if (le == 4)            //NUC OPTION :default is nuc_simple
 179                         {
 180                 read_matrix(nuc_simple);
 181                 gap_open = -4;
 182                 gap_ext = -0.25;
 183         }
 184
 185         else if (le == 160)  //PROT option: default is gonnet_160
 186                         {
 187                 if (TRACE)
 188                         printf("read matrix\n");
 189                 read_matrix(gonnet_160);
 190                 gap_open = -22;
 191                 gap_ext = -1;
 192
 193                 read_normalized_matrix(normalized_blosum_30); // add by YE Yongtao
 194         } else if (le == 1000) {  //Error handling
 195                 printf("Error: enter a valid matrix type\n");
 196                 exit(1);
 197                 //additional matrices can only be lower triangular
 198         }
 199
 200         //now override the gapopen and gapext
 201         if (argument.gapopen != 0.0 || argument.gapext != 0.00)
 202
 203         {
 204                 gap_open = -argument.gapopen;
 205                 gap_ext = -argument.gapext;
 206         }
 207
 208         if (TRACE)
 209                 printf("%f %f %f %d\n", argument.T, gap_open, gap_ext, le);
 210
 211         argument.gapopen = gap_open;
 212         argument.gapext = gap_ext;
 213         argument.opt = 'P';
 214
 215 }