00001 #include <ctype.h>
00002 #include "utils.h"
00003 #include "fold_vars.h"
00004
00005 #define NBASES 8
00006
00007
00008 static const char Law_and_Order[] = "_ACGUTXKI";
00009 static int BP_pair[NBASES][NBASES]=
00010
00011 {{ 0, 0, 0, 0, 0, 0, 0, 0},
00012 { 0, 0, 0, 0, 5, 0, 0, 5},
00013 { 0, 0, 0, 1, 0, 0, 0, 0},
00014 { 0, 0, 2, 0, 3, 0, 0, 0},
00015 { 0, 6, 0, 4, 0, 0, 0, 6},
00016 { 0, 0, 0, 0, 0, 0, 2, 0},
00017 { 0, 0, 0, 0, 0, 1, 0, 0},
00018 { 0, 6, 0, 0, 5, 0, 0, 0}};
00019
00020 #define MAXALPHA 20
00021
00022 static short alias[MAXALPHA+1];
00023 static int pair[MAXALPHA+1][MAXALPHA+1];
00024
00025 static int rtype[8] = {0, 2, 1, 4, 3, 6, 5, 7};
00026
00027 #ifdef _OPENMP
00028 #pragma omp threadprivate(Law_and_Order, BP_pair, alias, pair, rtype)
00029 #endif
00030
00031
00032 #define ENCODE(c) encode_char(c)
00033
00034 static int encode_char(char c) {
00035
00036 int code;
00037 if (energy_set>0) code = (int) (c-'A')+1;
00038 else {
00039 const char *pos;
00040 pos = strchr(Law_and_Order, c);
00041 if (pos==NULL) code=0;
00042 else code = (int) (pos-Law_and_Order);
00043 if (code>5) code = 0;
00044 if (code>4) code--;
00045 }
00046 return code;
00047 }
00048
00049
00050
00051 extern char *nonstandards;
00052 extern void nrerror(const char message[]);
00053 static void make_pair_matrix(void)
00054 {
00055 int i,j;
00056
00057 if (energy_set==0) {
00058 for (i=0; i<5; i++) alias[i] = (short) i;
00059 alias[5] = 3;
00060 alias[6] = 2;
00061 alias[7] = 0;
00062 for (i=0; i<NBASES; i++) {
00063 for (j=0; j<NBASES; j++)
00064 pair[i][j] = BP_pair[i][j];
00065 }
00066 if (noGU) pair[3][4] = pair[4][3] =0;
00067 if (nonstandards!=NULL) {
00068 for (i=0; i<(int)strlen(nonstandards); i+=2)
00069 pair[encode_char(nonstandards[i])]
00070 [encode_char(nonstandards[i+1])]=7;
00071 }
00072 for (i=0; i<NBASES; i++) {
00073 for (j=0; j<NBASES; j++)
00074 rtype[pair[i][j]] = pair[j][i];
00075 }
00076 } else {
00077 for (i=0; i<=MAXALPHA; i++) {
00078 for (j=0; j<=MAXALPHA; j++)
00079 pair[i][j] = 0;
00080 }
00081 if (energy_set==1) {
00082 for (i=1; i<MAXALPHA;) {
00083 alias[i++] = 3;
00084 alias[i++] = 2;
00085 }
00086 for (i=1; i<MAXALPHA; i++) {
00087 pair[i][i+1] = 2;
00088 i++;
00089 pair[i][i-1] = 1;
00090 }
00091 }
00092 else if (energy_set==2) {
00093 for (i=1; i<MAXALPHA;) {
00094 alias[i++] = 1;
00095 alias[i++] = 4;
00096 }
00097 for (i=1; i<MAXALPHA; i++) {
00098 pair[i][i+1] = 5;
00099 i++;
00100 pair[i][i-1] = 6;
00101 }
00102 }
00103 else if (energy_set==3) {
00104 for (i=1; i<MAXALPHA-2; ) {
00105 alias[i++] = 3;
00106 alias[i++] = 2;
00107 alias[i++] = 1;
00108 alias[i++] = 4;
00109 }
00110 for (i=1; i<MAXALPHA-2; i++) {
00111 pair[i][i+1] = 2;
00112 i++;
00113 pair[i][i-1] = 1;
00114 i++;
00115 pair[i][i+1] = 5;
00116 i++;
00117 pair[i][i-1] = 6;
00118 }
00119 }
00120 else nrerror("What energy_set are YOU using??");
00121 for (i=0; i<=MAXALPHA; i++) {
00122 for (j=0; j<=MAXALPHA; j++)
00123 rtype[pair[i][j]] = pair[j][i];
00124 }
00125 }
00126 }
00127
00128 static short *encode_sequence(const char *sequence, short how){
00129 unsigned int i,l = (unsigned int)strlen(sequence);
00130 short *S = (short *) space(sizeof(short)*(l+2));
00131
00132 switch(how){
00133
00134 case 0: for(i=1; i<=l; i++)
00135 S[i]= (short) encode_char(toupper(sequence[i-1]));
00136 S[l+1] = S[1];
00137 S[0] = (short) l;
00138 break;
00139
00140 case 1: for(i=1; i<=l; i++)
00141 S[i] = alias[(short) encode_char(toupper(sequence[i-1]))];
00142 S[l+1] = S[1];
00143 S[0] = S[l];
00144 break;
00145 }
00146
00147 return S;
00148 }