binaries/src/ViennaRNA/H/pair_mat.h

   1 #include <ctype.h>
   2 #include "utils.h"
   3 #include "fold_vars.h"
   4
   5 #define NBASES 8
   6 /*@notnull@*/
   7
   8 static const char Law_and_Order[] = "_ACGUTXKI";
   9 static int BP_pair[NBASES][NBASES]=
  10 /* _  A  C  G  U  X  K  I */
  11 {{ 0, 0, 0, 0, 0, 0, 0, 0},
  12  { 0, 0, 0, 0, 5, 0, 0, 5},
  13  { 0, 0, 0, 1, 0, 0, 0, 0},
  14  { 0, 0, 2, 0, 3, 0, 0, 0},
  15  { 0, 6, 0, 4, 0, 0, 0, 6},
  16  { 0, 0, 0, 0, 0, 0, 2, 0},
  17  { 0, 0, 0, 0, 0, 1, 0, 0},
  18  { 0, 6, 0, 0, 5, 0, 0, 0}};
  19
  20 #define MAXALPHA 20       /* maximal length of alphabet */
  21
  22 static short alias[MAXALPHA+1];
  23 static int pair[MAXALPHA+1][MAXALPHA+1];
  24 /* rtype[pair[i][j]]:=pair[j][i] */
  25 static int rtype[8] = {0, 2, 1, 4, 3, 6, 5, 7};
  26
  27 #ifdef _OPENMP
  28 #pragma omp threadprivate(Law_and_Order, BP_pair, alias, pair, rtype)
  29 #endif
  30
  31 /* for backward compatibility */
  32 #define ENCODE(c) encode_char(c)
  33
  34 static int encode_char(char c) {
  35   /* return numerical representation of base used e.g. in pair[][] */
  36   int code;
  37   if (energy_set>0) code = (int) (c-'A')+1;
  38   else {
  39     const char *pos;
  40     pos = strchr(Law_and_Order, c);
  41     if (pos==NULL) code=0;
  42     else code = (int) (pos-Law_and_Order);
  43     if (code>5) code = 0;
  44     if (code>4) code--; /* make T and U equivalent */
  45   }
  46   return code;
  47 }
  48
  49 /*@+boolint +charint@*/
  50 /*@null@*/
  51 extern char *nonstandards;
  52 extern void   nrerror(const char message[]);
  53 static void make_pair_matrix(void)
  54 {
  55    int i,j;
  56
  57    if (energy_set==0) {
  58       for (i=0; i<5; i++) alias[i] = (short) i;
  59       alias[5] = 3; /* X <-> G */
  60       alias[6] = 2; /* K <-> C */
  61       alias[7] = 0; /* I <-> default base '@' */
  62       for (i=0; i<NBASES; i++) {
  63           for (j=0; j<NBASES; j++)
  64             pair[i][j] = BP_pair[i][j];
  65       }
  66       if (noGU) pair[3][4] = pair[4][3] =0;
  67       if (nonstandards!=NULL) {  /* allow nonstandard bp's */
  68          for (i=0; i<(int)strlen(nonstandards); i+=2)
  69             pair[encode_char(nonstandards[i])]
  70               [encode_char(nonstandards[i+1])]=7;
  71       }
  72       for (i=0; i<NBASES; i++) {
  73           for (j=0; j<NBASES; j++)
  74            rtype[pair[i][j]] = pair[j][i];
  75       }
  76    } else {
  77       for (i=0; i<=MAXALPHA; i++) {
  78          for (j=0; j<=MAXALPHA; j++)
  79             pair[i][j] = 0;
  80       }
  81       if (energy_set==1) {
  82          for (i=1; i<MAXALPHA;) {
  83             alias[i++] = 3;  /* A <-> G */
  84             alias[i++] = 2;  /* B <-> C */
  85          }
  86          for (i=1; i<MAXALPHA; i++) {
  87             pair[i][i+1] = 2;    /* AB <-> GC */
  88             i++;
  89             pair[i][i-1] = 1;    /* BA <-> CG */
  90          }
  91       }
  92       else if (energy_set==2) {
  93         for (i=1; i<MAXALPHA;) {
  94             alias[i++] = 1;  /* A <-> A*/
  95             alias[i++] = 4;  /* B <-> U */
  96          }
  97          for (i=1; i<MAXALPHA; i++) {
  98             pair[i][i+1] = 5;    /* AB <-> AU */
  99             i++;
 100             pair[i][i-1] = 6;    /* BA <-> UA */
 101          }
 102       }
 103       else if (energy_set==3) {
 104         for (i=1; i<MAXALPHA-2; ) {
 105           alias[i++] = 3;  /* A <-> G */
 106           alias[i++] = 2;  /* B <-> C */
 107           alias[i++] = 1;  /* C <-> A */
 108           alias[i++] = 4;  /* D <-> U */
 109         }
 110         for (i=1; i<MAXALPHA-2; i++) {
 111           pair[i][i+1] = 2;    /* AB <-> GC */
 112           i++;
 113           pair[i][i-1] = 1;    /* BA <-> CG */
 114           i++;
 115           pair[i][i+1] = 5;    /* CD <-> AU */
 116           i++;
 117           pair[i][i-1] = 6;    /* DC <-> UA */
 118         }
 119       }
 120       else nrerror("What energy_set are YOU using??");
 121       for (i=0; i<=MAXALPHA; i++) {
 122         for (j=0; j<=MAXALPHA; j++)
 123           rtype[pair[i][j]] = pair[j][i];
 124       }
 125    }
 126 }
 127
 128 static short *encode_sequence(const char *sequence, short how){
 129   unsigned int i,l = (unsigned int)strlen(sequence);
 130   short         *S = (short *) space(sizeof(short)*(l+2));
 131
 132   switch(how){
 133     /* standard encoding as always used for S */
 134     case 0:   for(i=1; i<=l; i++) /* make numerical encoding of sequence */
 135                 S[i]= (short) encode_char(toupper(sequence[i-1]));
 136               S[l+1] = S[1];
 137               S[0] = (short) l;
 138               break;
 139     /* encoding for mismatches of nostandard bases (normally used for S1) */
 140     case 1:   for(i=1; i<=l; i++)
 141                 S[i] = alias[(short) encode_char(toupper(sequence[i-1]))];
 142               S[l+1] = S[1];
 143               S[0] = S[l];
 144               break;
 145   }
 146
 147   return S;
 148 }