sources/multicoil/PairCoilDiffer.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 PairCoilDifferDimension(int mode, Sequence sequence,
  16                              char lib[MAXFUNCTNUM], double *maxscore,
  17                              char preproc_offsets[MAXSEQLEN],
  18                              double preproc_scores[MAXSEQLEN][POSNUM+1],
  19                   double all_scores[MAX_NUM_SCORE_DIM][MAXSEQLEN][POSNUM+1],
  20                   char all_offsets[MAX_NUM_SCORE_DIM][MAXSEQLEN], int dimension,
  21                              int offset_to_use,
  22                              int number_differ_score_dim, double bound)
  23 {
  24   int i,j, pass,pos;
  25   double sc2scores[MAXSEQLEN][POSNUM];
  26   char diff_offsets[MAXSEQLEN];
  27   extern int which_differentiator;
  28   extern int differ_window_length;
  29   int number_passes=1;
  30   int start_coil =-1;
  31
  32   double coil_score[POSNUM];
  33
  34   /** which_differentiator = 1 means do window score starting at residue**/
  35   /** If it is 2 means find both max and min window.  ******/
  36   /** Passing 1 into scseqadj makes it return residue scores so the **/
  37   /** residue gets the window score for the window starting at that **/
  38   /** residue.  2 returns the max window score containing the residue **/
  39   /** and 3 returns the min window score. **/
  40
  41   WINDOW = differ_window_length;  /** TEMPORARY.  **/
  42
  43   if (which_differentiator ==2) number_passes=2;
  44   else number_passes=1;
  45
  46   for ( pass=0; pass<number_passes; pass++) {
  47     if (which_differentiator ==3 || which_differentiator==4) {
  48                                 /** get coiled regions from preproc. **/
  49       i=0;
  50       while(i<=sequence.seqlen) {
  51         if ( (start_coil >= 0) && ( (i==sequence.seqlen) ||
  52            (preproc_scores[i][7]<=bound) ||  /* End coil */
  53            ( (i>0) &&
  54             (preproc_offsets[i] != preproc_offsets[i-1])) )) { /*regshift*/
  55 /*        WINDOW= i-start_coil; */
  56 /*******
  57           scseqadj(lib,
  58                    &(sequence.seq[start_coil]), i-start_coil,
  59                    &(sc2scores[start_coil][0]),
  60                    &(diff_offsets[start_coil]),
  61                    maxscore, mode, 1);
  62 *******/
  63 /**          score_differentiator_window(lib,
  64                &(sequence.seq[start_coil]), i-start_coil,
  65                coil_score,mode, 1);
  66 **/
  67           score_differentiator_window(lib,
  68                &(sequence.seq[start_coil]), i-start_coil,
  69                coil_score,mode, 1);
  70
  71 /*       fprintf(stderr,"\n start_coil = %d, end_coil =%d",start_coil, i-1); */
  72
  73           for(j=start_coil; j<i; j++)
  74             for(pos=0;pos<POSNUM;pos++)
  75               if (which_differentiator==3)
  76                 sc2scores[j][pos]= coil_score[pos];
  77               else if (which_differentiator==4)
  78                 sc2scores[j][pos] = coil_score[pos]/(i-start_coil);
  79
  80           start_coil = -1;
  81
  82         }
  83
  84         if ( (preproc_scores[i][7] > bound) && (start_coil ==-1))
  85           start_coil =i;
  86         i++;
  87       }
  88     }
  89     else
  90       scseqadj(lib,
  91                sequence.seq, sequence.seqlen,
  92                sc2scores,
  93                diff_offsets,
  94                maxscore, mode, pass + which_differentiator);
  95                                /* Get differentiator score.*/
  96                               /* pass + which_diff = 2 is max_window_score **/
  97                               /*                     3 is min_window_score **/
  98
  99     for (i=0; i<sequence.seqlen; ++i) {
 100       /* Use the offset given by the preproc if offset_to_use is -1 or 7 */
 101       if ((offset_to_use == -1) || (offset_to_use ==7))  {
 102         all_scores[dimension+ pass*number_differ_score_dim/2][i][7] =
 103           sc2scores[i][preproc_offsets[i]];
 104         all_offsets[dimension+ pass*number_differ_score_dim/2][i] =
 105           preproc_offsets[i];
 106       }
 107       else {    /*  Offset to use is 0-6. */
 108         all_scores[dimension+ pass*number_differ_score_dim/2][i][7] =
 109           sc2scores[i][offset_to_use];
 110         all_offsets[dimension+ pass*number_differ_score_dim/2][i]=
 111           offset_to_use;
 112       }
 113
 114
 115       /* do other registers */
 116       for (j=0; j<7; ++j)  { /** sc2scores[i][j] is the score for position
 117                              i+j%POSNUM  **/
 118         all_scores[dimension+ pass*number_differ_score_dim/2][i][(i+j)%POSNUM]=
 119           sc2scores[i][j];
 120       }
 121     }
 122   }
 123
 124 }
 125
 126
 127
 128
 129 /******************************************************************/
 130 void zero_out_non_coils(
 131                      double all_scores[MAX_NUM_SCORE_DIM][MAXSEQLEN][POSNUM+1],
 132                      double preproc_like[MAXSEQLEN][POSNUM+1],
 133                      int seqlen, double bound)
 134 {
 135
 136   int res, reg, olig;
 137
 138   for(olig=0; olig<2; olig++)
 139     for (res=0; res<seqlen; res++)
 140       for (reg=0;reg<=POSNUM;reg++)
 141         if (preproc_like[res][7] > bound)
 142           all_scores[olig][res][reg] = all_scores[olig+3][res][reg];
 143         else
 144           all_scores[olig][res][reg] =0;
 145 }
 146
 147
 148 void zero_out_bad_windows(int differ_window_length,
 149                      double all_scores[MAX_NUM_SCORE_DIM][MAXSEQLEN][POSNUM+1],
 150                      double preproc_like[MAXSEQLEN][POSNUM+1],
 151                      int seqlen, double bound)
 152 {
 153   int i,j;
 154   int last_bad_window_start=-1;
 155   int olig_state, pos;
 156 /*** For now just look at score of first and last residue in window to **/
 157 /**  see if they are above bound.  To really do it right should look at */
 158 /** min score in window (using min_over_windows on preproc_like), but   */
 159 /** the format expected by min_over_windows would require the POSNUM+1  */
 160 /** indice of preproc_like be in the first, not second position....     */
 161
 162   for (i=0; i<seqlen; i++)
 163     if (preproc_like[i][7] < bound)  { /** Zero out windows. **/
 164       for (j=MAX(last_bad_window_start+1,i-differ_window_length); j<=i; j++)
 165         for(olig_state=0; olig_state<2; olig_state++)
 166           for(pos=0; pos<=POSNUM; pos++)
 167             all_scores[olig_state][j][pos]=0;
 168
 169       last_bad_window_start=i;
 170     }
 171     else   /* Do this in case bound has changed....using stored value in +3 */
 172       for(olig_state=0; olig_state<2; olig_state++)
 173         for(pos=0; pos<=POSNUM; pos++)
 174           all_scores[olig_state][i][pos]=all_scores[olig_state+3][i][pos];
 175 }
 176
 177
 178
 179