sources/multicoil/PairCoilDiffer_cuts.c

   1 #include <stdio.h>
   2 #include <stdlib.h>
   3 #include <math.h>
   4 #include <string.h>
   5 #include "scio.h"
   6 #include "switches.h"
   7 #include "likelihood.h"
   8 #include "sc.h"
   9 #include "scconst.h"
  10 #include "options.h"
  11 #include "compute_like.h"
  12 #include "interface.h"
  13 #include "stats.h"
  14
  15 void average_score_over_coils2(Sequence seq,
  16                                double scores[MAXSEQLEN][POSNUM+1],
  17                               double PreprocLike[MAXSEQLEN][POSNUM+1],
  18                               int CoilDiffMethod, int offset_to_use,
  19                               char offsets[MAXSEQLEN],
  20                               int start_at_reg_shift)
  21 {
  22   int i=0, j,k;
  23   int start_coil[POSNUM+1];
  24   double total_coil_score[POSNUM+1];
  25   double avg_coil_score[POSNUM+1];
  26   int regist;
  27
  28   double total_coil_length[POSNUM +1];  /* A weighted length. */
  29
  30   for (j=0; j<POSNUM+1; j++) {
  31     start_coil[j] = -1;
  32     total_coil_score[j]=0;
  33     avg_coil_score[j]=0;
  34     total_coil_length[j]=0;
  35   }
  36
  37   if (offset_to_use == -1 ) offset_to_use = 7;  /* Adjust "all" offset */
  38                                               /* to max offset locally. */
  39
  40   while (i<=seq.seqlen) {
  41     for (k=0; k<=POSNUM; k++) {
  42       if (k!=7) regist = (i+k)%POSNUM;
  43       else regist =7;
  44
  45       /*  Consider < 1/10 percent to be zero */
  46       if ( ( PreprocLike[i][regist] < .001) || (i==seq.seqlen) ||
  47           ( (i>0) && (offsets[i] != offsets[i-1])  && start_at_reg_shift &&
  48            (PreprocLike[i-1][7]>= .001) ) )  {
  49       /** If not in a coil (value of -HUGE_VAL signals not in coil).  */
  50         if (start_coil[k] != -1)  {  /* Ended coil at i-1. */
  51           avg_coil_score[k] = total_coil_score[k]/total_coil_length[k];
  52           for  (j=start_coil[k]; j<i; j++) {
  53             if (k!=7)
  54               scores[j][(j+k)%POSNUM]= avg_coil_score[k];
  55             else
  56               scores[j][7]= avg_coil_score[k];
  57
  58           }
  59         }
  60         total_coil_score[k] =0;
  61 void average_score_over_coils2(Sequence seq,
  62                                double scores[MAXSEQLEN][POSNUM+1],
  63                               double PreprocLike[MAXSEQLEN][POSNUM+1],
  64                               int CoilDiffMethod, int offset_to_use,
  65                               char offsets[MAXSEQLEN],
  66                               int start_at_reg_shift)
  67 {
  68   int i=0, j,k;
  69   int start_coil[POSNUM+1];
  70   double total_coil_score[POSNUM+1];
  71   double avg_coil_score[POSNUM+1];
  72   int regist;
  73
  74   double total_coil_length[POSNUM +1];  /* A weighted length. */
  75
  76   for (j=0; j<POSNUM+1; j++) {
  77     start_coil[j] = -1;
  78     total_coil_score[j]=0;
  79     avg_coil_score[j]=0;
  80     total_coil_length[j]=0;
  81   }
  82
  83   if (offset_to_use == -1 ) offset_to_use = 7;  /* Adjust "all" offset */
  84                                               /* to max offset locally. */
  85
  86   while (i<=seq.seqlen) {
  87     for (k=0; k<=POSNUM; k++) {
  88       if (k!=7) regist = (i+k)%POSNUM;
  89       else regist =7;
  90
  91       /*  Consider < 1/10 percent to be zero */
  92       if ( ( PreprocLike[i][regist] < .001) || (i==seq.seqlen) ||
  93           ( (i>0) && (offsets[i] != offsets[i-1])  && start_at_reg_shift &&
  94            (PreprocLike[i-1][7]>= .001) ) )  {
  95       /** If not in a coil (value of -HUGE_VAL signals not in coil).  */
  96         if (start_coil[k] != -1)  {  /* Ended coil at i-1. */
  97           avg_coil_score[k] = total_coil_score[k]/total_coil_length[k];
  98           for  (j=start_coil[k]; j<i; j++) {
  99             if (k!=7)
 100               scores[j][(j+k)%POSNUM]= avg_coil_score[k];
 101             else
 102               scores[j][7]= avg_coil_score[k];
 103
 104           }
 105         }
 106         total_coil_score[k] =0;
 107         if (i != seq.seqlen)   /* Consider < 1/10 percent to be 0 */
 108           if  ( PreprocLike[i][regist] < .001)  {  /* Not in coil anymore. */
 109             scores[i][regist]=0;
 110             start_coil[k] = -1;
 111             total_coil_length[k] =0;
 112           }
 113           else {  /* Changing registers. */
 114             start_coil[k]= i;
 115             if (CoilDiffMethod==0) {
 116               total_coil_length[k] = 1;
 117               total_coil_score[k] = scores[i][regist];
 118             }
 119             else if (CoilDiffMethod ==1) {
 120               total_coil_length[k]= PreprocLike[i][regist];
 121               total_coil_score[k] = scores[i][regist]*
 122                 PreprocLike[i][regist];
 123             }
 124           }
 125       }   /** Ends if ending a coil.  **/
 126
 127       else {    /* If in a coil.  */
 128         if (CoilDiffMethod ==0) {   /* Average scores over coil residues. */
 129           total_coil_score[k] += scores[i][regist];
 130
 131           total_coil_length[k] += 1;
 132         }
 133         else if (CoilDiffMethod == 1) {  /* Weighted average scores.  */
 134           total_coil_score[k] += scores[i][regist]*
 135                  PreprocLike[i][regist];
 136           total_coil_length[k] += PreprocLike[i][regist];
 137         }
 138         if (start_coil[k] == -1)
 139           start_coil[k] =i;
 140       }
 141     }
 142
 143     i++;
 144   }
 145
 146   /*  Need the following hack to get the correct register in max off method.*/
 147   if ( (offset_to_use == 7) && (!start_at_reg_shift))
 148     for (i=0; i<seq.seqlen; i++) {
 149       for (j=0; j<POSNUM; j++)
 150         scores[i][j]= -HUGE_VAL;
 151       scores[i][(i+offsets[i]) %POSNUM] = scores[i][7];
 152     }
 153 }
 154
 155
 156
 157 /*******************************************************************/
 158
 159 void PairCoilDiffer(int mode, Sequence sequence, char lib[MAXFUNCTNUM],
 160                       double *maxscore,char offsets[MAXSEQLEN],
 161                       double preproc_like[MAXSEQLEN][POSNUM+1],
 162                       int offset_to_use, int table,
 163                       double scores[MAXSEQLEN][POSNUM+1],
 164                       int Avg_Coil_Score, int CoilDiffMethod,
 165                       int start_coil_at_reg_shift)
 166
 167 {
 168   int i,j,k;
 169   double sc2scores[MAXSEQLEN][POSNUM];
 170   char diff_offsets[MAXSEQLEN];
 171   double total_coilscore[POSNUM+1];
 172   int coil_length[POSNUM +1];
 173   int like2or3;
 174
 175   extern long pfreqs[MAX_TABLE_NUMBER][AANUM][POSNUM];
 176   extern long ptotals[MAX_TABLE_NUMBER][POSNUM];
 177   extern long pfreqp[MAX_TABLE_NUMBER][AANUM][AANUM][POSNUM][POSNUM];
 178   extern long ptotalp[MAX_TABLE_NUMBER][POSNUM][POSNUM];
 179
 180
 181   int ProlineFreeWindow= mode & PROLINE_FREE_MODE;
 182
 183   if (table ==0) like2or3= 2;
 184   else like2or3=3;
 185
 186   for (i=0; i<=POSNUM; i++) {
 187     coil_length[i] = 0;
 188     total_coilscore[i] =0;
 189   }
 190
 191
 192
 193   /* score the sequence using sc2seq */
 194   scseqadj(lib,
 195            sequence.seq, sequence.seqlen,
 196            sc2scores,
 197            diff_offsets,
 198            maxscore, mode, 1);  /* Get differentiator score.*/
 199
 200
 201
 202
 203
 204 /*  *maxscore= exp(*maxscore);   */
 205
 206   *maxscore =0;
 207
 208   for (i=0; i<sequence.seqlen; ++i) {
 209     /*  Do the most probabable register if offset_to_use is -1 or 7. */
 210     if ((offset_to_use == -1) || (offset_to_use ==7))
 211       if (like2or3 ==2)
 212         scores[i][7] = 1 - exp(sc2scores[i][offsets[i]]);
 213       else
 214         scores[i][7] =  exp(sc2scores[i][offsets[i]]);
 215
 216     else    /*  Offset to use is 0-6. */
 217       if (like2or3 ==  2)
 218         scores[i][7] = 1 - exp(sc2scores[i][offset_to_use]);
 219       else
 220         scores[i][7] = exp(sc2scores[i][offset_to_use]);
 221                 /*****Trimer like is 1 - dimerlike ***/
 222
 223     if (preproc_like[i][7] <.001)   /* Consider < 1/10 percent to be 0 */
 224       scores[i][7]=0;  /* Don't score non-coiled regions. */
 225
 226
 227
 228     for (j=0; j<7; ++j)  { /** sc2scores[i][j] is the score for position
 229                              i+j%POSNUM  **/
 230       if (preproc_like[i][(i+j)%POSNUM]> 0) {
 231         if (like2or3==2)
 232           scores[i][(i+j)%POSNUM] = 1 - exp(sc2scores[i][j]);
 233         else
 234           scores[i][(i+j)%POSNUM] = exp(sc2scores[i][j]);
 235       }
 236       else
 237         scores[i][(i+j)%POSNUM]= 0;
 238     }
 239
 240   }
 241
 242
 243   if (Avg_Coil_Score)
 244     average_score_over_coils2(sequence, scores, preproc_like,
 245                               CoilDiffMethod, offset_to_use,offsets,
 246                               start_coil_at_reg_shift);
 247
 248
 249
 250   for (i=0; i<sequence.seqlen; i++)
 251     if (scores[i][7] > *maxscore) *maxscore = scores[i][7];
 252
 253
 254 }
 255
 256
 257