00001 #include <math.h>
00002 #include <stdio.h>
00003 #include <stdlib.h>
00004 #include <string.h>
00005 #include <ctype.h>
00006 #include <assert.h>
00007
00008 #include "open_prospect.h"
00009
00010
00012
00014
00015
00016 void append_to_str(char **str1, char *str2) {
00017 int str_len1, str_len2,i;
00018 char *tmp;
00019
00020 if (*str1 == NULL)
00021 str_len1 = 0;
00022 else
00023 str_len1 = strlen(*str1);
00024
00025 if (str2 == NULL)
00026 str_len2 = 0;
00027 else
00028 str_len2 = strlen(str2);
00029
00030 tmp = (char *)malloc(sizeof(char) * (str_len1 + str_len2 + 1));
00031 assert( tmp != NULL);
00032 for (i = 0; i < str_len1; i++) {
00033 tmp[i] = (*str1)[i];
00034 }
00035 for (i = 0; i < str_len2; i++) {
00036 tmp[i+str_len1] = str2[i];
00037 }
00038
00039 tmp[str_len1 + str_len2] = 0;
00040 if (*str1 != NULL)
00041 free(*str1);
00042 *str1 = tmp;
00043 }
00044
00045
00046
00047 AlignmentStruct::AlignmentStruct(const AlignmentStruct &rhs) {
00048 method[0] = 0;
00049 Copy(rhs);
00050 }
00051
00052 AlignmentStruct & AlignmentStruct::operator=(const AlignmentStruct &rhs) {
00053 Free();
00054 Copy(rhs);
00055 return *this;
00056 }
00057
00058
00059
00060 void AlignmentStruct::Copy(const AlignmentStruct &rhs) {
00061 method[0] = 0;
00062 strncpy( method, rhs.method, sizeof( rhs.method ) );
00063 target_len = rhs.target_len;
00064 template_len = rhs.template_len;
00065 target_align = (int *)malloc(sizeof(int) * target_len);
00066 template_align = (int *)malloc(sizeof(int) * template_len);
00067 memcpy( target_align, rhs.target_align, sizeof(int) * target_len );
00068 memcpy( template_align, rhs.template_align, sizeof(int) * template_len );
00069 frag_count = 0;
00070 frag_start = NULL;
00071 frag_stop = NULL;
00072 }
00073
00074
00075 AlignmentStruct::AlignmentStruct() {
00076 target_align = NULL;
00077 template_align = NULL;
00078 frag_start = NULL;
00079 frag_stop = NULL;
00080 method[0] = 0;
00081 frag_count = 0;
00082 }
00083
00084 AlignmentStruct::~AlignmentStruct() {
00085 Free();
00086 }
00087
00088 void AlignmentStruct::Free(void) {
00089 if (target_align != NULL)
00090 free(target_align);
00091 target_align = NULL;
00092 if (template_align != NULL)
00093 free(template_align);
00094 template_align = NULL;
00095 if (frag_start != NULL)
00096 free(frag_start);
00097 frag_start = NULL;
00098 if (frag_stop != NULL)
00099 free(frag_stop);
00100 frag_stop = NULL;
00101 }
00102
00103
00104 int AlignmentStruct::GetFragCount(void) {
00105 if ( target_align != NULL && frag_start == NULL)
00106 CalcFrag();
00107 return frag_count;
00108 }
00109
00110
00111 int AlignmentStruct::GetFragStart(int i) {
00112 if ( target_align != NULL && frag_start == NULL)
00113 CalcFrag();
00114 if ( i >= 0 && i < frag_count )
00115 return frag_start[i];
00116 return -1;
00117 }
00118
00119 int AlignmentStruct::GetFragStop(int i) {
00120 if ( target_align != NULL && frag_stop == NULL)
00121 CalcFrag();
00122 if ( i >= 0 && i < frag_count )
00123 return frag_stop[i];
00124 return -1;
00125 }
00126
00127
00128 void AlignmentStruct::CalcFrag(void) {
00129 if ( target_align == NULL)
00130 return;
00131 frag_count = 0;
00132
00133 #if 0 //trying out different method of frag counting, trying to figure out one that works all the time
00134 for (int cycle = 0; cycle < 2; cycle++) {
00135 int cur_frag = 0;
00136
00137 int last_pos = -1;
00138 int start_pos = -100, end_pos = -100;
00139 for (int j = 0; j < target_len; j++) {
00140 if ( (target_align[j] == -1 || j == target_len-1) && last_pos != -1) {
00141
00142 if ( target_len-1 == j && target_align[j] != -1 )
00143 end_pos = target_len-1;
00144 else
00145 end_pos = j-1;
00146 if (cycle == 0) {
00147 frag_count++;
00148 } else {
00149 frag_start[cur_frag] = start_pos;
00150 frag_stop[cur_frag] = end_pos + 1;
00151 cur_frag++;
00152 }
00153 }
00154 if (target_align[j] != -1 && (target_align[j] != last_pos + 1 || target_align[j] == 0) ) {
00155 if (last_pos != -1) {
00156 end_pos = j-1;
00157 if ( cycle == 0 ) {
00158 frag_count++;
00159 } else {
00160 frag_start[cur_frag] = start_pos;
00161 frag_stop[cur_frag] = end_pos + 1;
00162 cur_frag++;
00163 }
00164 }
00165 start_pos = j;
00166 last_pos = target_align[j];
00167 }
00168 last_pos = target_align[j];
00169 }
00170 if (cycle == 0) {
00171 frag_start = (int *)malloc(sizeof(int) * frag_count);
00172 frag_stop = (int *)malloc(sizeof(int) * frag_count);
00173 }
00174 }
00175 #else
00176 int last_pos = -10;
00177 for (int i = 0 ;i < target_len; i++) {
00178 if ( target_align[i] != -1 ) {
00179 if ( target_align[i] != last_pos + 1 )
00180 frag_count++;
00181 last_pos = target_align[i];
00182 } else {
00183 last_pos = -10;
00184 }
00185 }
00186 frag_start = (int *)malloc(sizeof(int) * frag_count);
00187 frag_stop = (int *)malloc(sizeof(int) * frag_count);
00188 int cur_frag = 0;
00189 for (int i = 0 ; i < target_len; ) {
00190 if ( target_align[i] != -1 ) {
00191 int j = i+1;
00192 while ( j < target_len && target_align[j] == target_align[j-1]+1)
00193 j++;
00194 frag_start[ cur_frag ] = i;
00195 frag_stop [ cur_frag ] = j;
00196 i = j;
00197 cur_frag++;
00198 } else {
00199 i++;
00200 }
00201 }
00202 assert( cur_frag == frag_count );
00203 #endif
00204 }
00205
00206
00207
00208
00209 void AlignmentStruct::GetEnergies( ProspectParam *param,
00210 TemplateStruct *template_data, TargetStruct *target_data,
00211 EnergyArray energy, long mask) {
00212 int gap_open;
00213 int target_start = 0;
00214 int target_stop = target_data->len;
00215
00216
00217 for (int i = 0; i < e_count; i++)
00218 energy[i] = 0;
00219
00220 if ( target_data->twobody == NULL )
00221 target_data->Setup( param );
00222
00223
00224 if (! param->no_loop) {
00225 for (int i = target_start; i < target_stop; i++) {
00226 if (this->target_align[i] != -1) {
00227 if (this->target_align[i] < template_data->len) {
00228 if ( e_mask( e_mutation ) & mask )
00229 energy[e_mutation] += calc_mutation(param, template_data, target_data, this->target_align[i],i);
00230
00231 if ( e_mask( e_mutationlog ) & mask )
00232 energy[e_mutationlog] += calc_mutationlog(param, template_data, target_data, this->target_align[i],i);
00233
00234 if ( e_mask( e_mutationfull ) & mask )
00235 energy[e_mutationfull] += calc_mutationfull(param, template_data, target_data, this->target_align[i],i);
00236
00237 if ( e_mask( e_singleton ) & mask )
00238 energy[e_singleton] += calc_singleton(param, template_data, target_data,this->target_align[i],i);
00239
00240 if ( e_mask( e_sec_struct ) & mask )
00241 energy[e_sec_struct] += calc_sec_struct(param, template_data, target_data,this->target_align[i],i);
00242
00243 if ( e_mask( e_mutscore ) & mask )
00244 energy[e_mutscore] += calc_mutscore(param, template_data, target_data,this->target_align[i],i);
00245
00246 if ( e_mask( e_twobody_frozen ) & mask )
00247 energy[e_twobody_frozen] += calc_twobody_frozen(param, template_data, target_data,this->target_align[i],i);
00248
00249 if ( e_mask( e_contact ) & mask )
00250 energy[e_contact] += calc_contact_frozen(param, template_data, target_data,this->target_align[i],i);
00251 }
00252 for (int j = 0; j < e_count; j++)
00253 assert( !isnan(energy[j]));
00254 }
00255 }
00256 } else {
00257
00258 for (int i = 0; i < template_data->core_count; i++) {
00259 for (int j = 0; j < template_data->cores[i].len; j++) {
00260 int k = template_data->cores[i].head + j;
00261 int target_k = this->template_align[k];
00262 if (target_k >= target_start && target_k < target_stop && k != -1) {
00263 if ( e_mask( e_mutation ) & mask )
00264 energy[e_mutation] += calc_mutation(param, template_data, target_data, k, target_k);
00265
00266 if ( e_mask( e_mutationlog ) & mask )
00267 energy[e_mutationlog] += calc_mutationlog(param, template_data, target_data, k, target_k);
00268
00269 if ( e_mask( e_mutationfull ) & mask )
00270 energy[e_mutationfull] += calc_mutationfull(param, template_data, target_data, k, target_k);
00271
00272 if ( e_mask( e_singleton ) & mask )
00273 energy[e_singleton] += calc_singleton(param, template_data, target_data, k, target_k);
00274
00275 if ( e_mask( e_sec_struct ) & mask )
00276 energy[e_sec_struct] += calc_sec_struct(param, template_data, target_data, k, target_k);
00277
00278 if ( e_mask( e_mutscore ) & mask )
00279 energy[e_mutscore] += calc_mutscore(param, template_data, target_data, k, target_k);
00280
00281 if ( e_mask( e_twobody_frozen ) & mask )
00282 energy[e_twobody_frozen] += calc_twobody_frozen(param, template_data, target_data, k, target_k);
00283
00284 if ( e_mask( e_contact ) & mask )
00285 energy[e_contact] += calc_contact_frozen(param, template_data, target_data, k,i);
00286 }
00287 }
00288 }
00289 }
00290
00291
00292 if ( !param->no_loop ) {
00293 if ( e_mask( e_gapopen ) & mask ) {
00294 gap_open = 0;
00295 for (int i = 0; i < target_data->len; i++) {
00296 if ( this->target_align[i] == -1) {
00297 energy[e_gapconst] += 1;
00298 if (!gap_open) {
00299 energy[e_gapopen] += 1;
00300 gap_open = 1;
00301 }
00302 } else
00303 gap_open = 0;
00304 }
00305 gap_open = 0;
00306 for (int i = 0; i < template_data->len; i++) {
00307 if ( this->template_align[i] == -1) {
00308 energy[e_gapconst] += 1;
00309 if (!gap_open) {
00310 energy[e_gapopen] += 1;
00311 gap_open = 1;
00312 }
00313 } else
00314 gap_open = 0;
00315 }
00316 }
00317 }
00318
00319
00320 if (template_data->cores) {
00321 if ( e_mask( e_coredel ) & mask ) {
00322 for (int i = 0; i < template_data->core_count; i++) {
00323 for (int j = 0; j < template_data->cores[i].len; j++) {
00324 int k = template_data->cores[i].head + j;
00325 int target_k = this->template_align[k];
00326 if ( target_k == -1 ) {
00327 energy[ e_coredel ] += 1 / (float) template_data->cores[i].len;
00328 }
00329 }
00330 }
00331 }
00332 }
00333
00334 {
00335 int last_target_pos = 0;
00336 for (int i = 0; i < template_data->len; i++) {
00337 int target_k = this->template_align[i];
00338 if ( target_k == -1 ) {
00339 if ( e_mask( e_tsingledel ) & mask )
00340 energy[ e_tsingledel ] += calc_singledel( param, template_data, i );
00341 if ( e_mask( e_contactdel ) & mask )
00342 energy[ e_contactdel ] += calc_contactdel( param, template_data, i );
00343 if ( e_mask( e_qsingleins ) & mask )
00344 energy[ e_qsingleins ] += target_data->residue[ last_target_pos ].single_ins;
00345 } else {
00346 last_target_pos = target_k;
00347 }
00348 }
00349 }
00350
00351 {
00352
00353 int last_template_pos = 0;
00354 for ( int i = 0 ; i < target_data->len; i++ ) {
00355 int template_k = this->target_align[i];
00356 if ( template_k == -1 ) {
00357 if ( e_mask( e_tsingleins ) & mask )
00358 energy[ e_tsingleins ] += calc_singleins( param, template_data, last_template_pos );
00359 if ( e_mask( e_qsingledel ) & mask )
00360 energy[ e_qsingledel ] += target_data->residue[i].single_del;
00361 } else {
00362 last_template_pos = template_k;
00363 }
00364 }
00365 }
00366
00367
00368 int template_res_num_1,template_res_num_2;
00369 int target_res_num_1,target_res_num_2;
00370
00371
00372 if ( param->full_twobody ) {
00373 for (int i = 0; i < template_data->len; i++) {
00374 target_res_num_1 = this->template_align[i];
00375 if ( target_res_num_1 != -1 ) {
00376 for (int j = 0; j < template_data->len; j++) {
00377 target_res_num_2 = this->template_align[i];
00378 if ( target_res_num_2 != -1 ) {
00379 if ( e_mask( e_twobody ) & mask )
00380 energy[e_twobody] += calc_twobody(param,template_data,target_data,i,j,
00381 target_res_num_1,target_res_num_2);
00382 if ( e_mask( e_dfire ) & mask )
00383 energy[e_dfire] += calc_dfire(param,template_data,target_data,i,j,
00384 target_res_num_1,target_res_num_2);
00385 }
00386 }
00387 }
00388 }
00389 } else {
00390 for (int i = 0; i < template_data->core_count; i++) {
00391 for (int j = i+1; j < template_data->core_count; j++) {
00392 char found = 0;
00393 if ( j != i+1 ) {
00394 for (int k=0; k < template_data->twobody_count && !found; k++) {
00395 if ( (template_data->twobody_interactions[k].core_number[0] == i &&
00396 template_data->twobody_interactions[k].core_number[1] == j) ||
00397 (template_data->twobody_interactions[k].core_number[0] == j &&
00398 template_data->twobody_interactions[k].core_number[1] == i) ) {
00399 found = 1;
00400 }
00401 }
00402 } else {
00403 found = 1;
00404 }
00405 if ( found ) {
00406 pValType twobody = 0, dfire = 0;
00407 for (int x = 0; x < template_data->cores[i].len; x++) {
00408 for (int y = 0; y < template_data->cores[j].len; y++) {
00409
00410
00411 template_res_num_1 = template_data->cores[i].head + x;
00412 template_res_num_2 = template_data->cores[j].head + y;
00413 target_res_num_1 = this->template_align[template_res_num_1];
00414 target_res_num_2 = this->template_align[template_res_num_2];
00415 if (target_res_num_1 >= target_start && target_res_num_2 >= target_start &&
00416 target_res_num_1 < target_stop && target_res_num_2 < target_stop) {
00417 if ( e_mask( e_twobody ) & mask )
00418 twobody += calc_twobody(param,template_data,target_data,template_res_num_1,template_res_num_2,
00419 target_res_num_1,target_res_num_2);
00420 if ( e_mask( e_dfire ) & mask )
00421 dfire += calc_dfire(param,template_data,target_data,template_res_num_1,template_res_num_2,
00422 target_res_num_1,target_res_num_2);
00423
00424 }
00425 }
00426 }
00427 energy[e_twobody] += twobody;
00428 energy[e_dfire] += dfire;
00429 }
00430 }
00431 }
00432 }
00433
00434
00435 if ( e_mask( e_distconst ) & mask ) {
00436 for (int i = 0; i < target_data->distconst_count; i++) {
00437 target_res_num_1 = target_data->distconst[i].aa[0];
00438 target_res_num_2 = target_data->distconst[i].aa[1];
00439 if (target_res_num_1 >= target_start && target_res_num_2 >= target_start &&
00440 target_res_num_1 < target_stop && target_res_num_2 < target_stop) {
00441 template_res_num_1 = this->target_align[ target_res_num_1 ];
00442 template_res_num_2 = this->target_align[ target_res_num_2 ];
00443 if (template_res_num_1 != -1 && template_res_num_2 != -1)
00444 energy[e_distconst]+= calc_distconst(param,template_data,target_data,template_res_num_1,template_res_num_2,
00445 target_res_num_1,target_res_num_2);
00446 }
00447 }
00448 }
00449 }
00450
00451
00452
00453
00454
00455
00456
00457 char * AlignmentStruct::Write_XML(ProspectParam *param,
00458 WeightArray weight,
00459 TemplateStruct *template_data,
00460 TargetStruct *target_data,
00461 ScoreStruct *score) {
00462
00463 char *outxml = NULL, buffer[4048], *template_str, *target_str,
00464 *template_ss_str, *target_ss_str, *align_str, *core_str;
00465 long i;
00466
00467 sprintf(buffer, "<threading template=\"%s\" method='%s'>\n", template_data->template_name, this->method);
00468 append_to_str(&outxml, buffer);
00469
00470 sprintf(buffer,
00471 "\t<templateInfo>\n<pdbcode>%s</pdbcode>\n<pdbchain>%c</pdbchain></templateInfo>\n",
00472 template_data->pdbcode, template_data->pdbchain);
00473 append_to_str(&outxml, buffer);
00474
00475 EnergyArray energy;
00476
00477 GetEnergies( param, template_data, target_data, energy);
00478
00479 if ( param->no_loop )
00480 sprintf(buffer, "\t<energy no_loop='yes'>\n");
00481 else
00482 sprintf(buffer, "\t<energy>\n");
00483 append_to_str(&outxml, buffer);
00484
00485 for (i = 0; i < e_count; i++) {
00486 if ( energy[i] != 0.0 ) {
00487 sprintf(buffer, "\t\t<%s e=\"%g\">%g</%s>\n", e_names[i], energy[i], weight[i] * energy[i], e_names[i]);
00488 append_to_str(&outxml, buffer);
00489 }
00490 }
00491
00492 sprintf(buffer, "\t</energy>\n");
00493
00494 append_to_str(&outxml, buffer);
00495
00496 if ( score && score->mask ) {
00497 sprintf(buffer, "\t<score>\n");
00498 append_to_str(&outxml, buffer);
00499 for (i = 0; i < ScoreStruct::SCORE_COUNT; i++) {
00500 if ( score->mask & (0x01<<i) ) {
00501 if ( i == ScoreStruct::SCORE_Z || i == ScoreStruct::SCORE_ZFULL ) {
00502 sprintf(buffer, "\t\t<%s mean=\"%g\" sd=\"%g\">%g</%s>\n",
00503 ScoreStruct::GetScoreName(i),
00504 (float)score->z_mean[i],
00505 (float)score->z_sd[i],
00506 (float)score->values[i],
00507 ScoreStruct::GetScoreName(i) );
00508
00509 } else {
00510 sprintf(buffer, "\t\t<%s>%g</%s>\n",
00511 ScoreStruct::GetScoreName(i),
00512 (float)score->values[i],
00513 ScoreStruct::GetScoreName(i));
00514 }
00515 append_to_str(&outxml, buffer);
00516 }
00517 }
00518 if ( score->mask & ScoreStruct::DO_SCORE_ZENERGY ) {
00519 for ( i = ScoreStruct::SCORE_ZENERGY_START; i < ScoreStruct::SCORE_ZENERGY_END; i++) {
00520 sprintf(buffer, "\t\t<%s mean=\"%g\" sd=\"%g\">%g</%s>\n",
00521 ScoreStruct::GetScoreName(i),
00522 (float)score->z_mean[i],
00523 (float)score->z_sd[i],
00524 (float)score->values[i],
00525 ScoreStruct::GetScoreName(i) );
00526 append_to_str(&outxml, buffer);
00527 }
00528 }
00529 sprintf(buffer, "\t</score>\n");
00530 append_to_str(&outxml, buffer);
00531 }
00532
00533 sprintf(buffer, "\t<weights>\n");
00534 append_to_str(&outxml, buffer);
00535 for (i = 0; i < e_count; i++) {
00536 if ( weight[i] != 0.0 ) {
00537 sprintf(buffer, "\t\t<%s>%g</%s>\n", e_names[i], weight[i], e_names[i]);
00538 append_to_str(&outxml, buffer);
00539 }
00540 }
00541 sprintf(buffer, "\t</weights>\n");
00542 append_to_str(&outxml, buffer);
00543 #if 0
00544 int alignment_strlen = target_data->len + template_data->len;
00545
00546
00547
00548
00549
00550
00551
00552
00553
00554
00555
00556
00557
00558
00559
00560 template_str = (char *)malloc(sizeof(char) * (alignment_strlen+1) );
00561 template_ss_str = (char *)malloc(sizeof(char) * (alignment_strlen+1) );
00562 target_str = (char *)malloc(sizeof(char) * (alignment_strlen+1) );
00563 target_ss_str = (char *)malloc(sizeof(char) * (alignment_strlen+1) );
00564 align_str = (char *)malloc(sizeof(char) * (alignment_strlen+1) );
00565 core_str = (char *)malloc(sizeof(char) * (alignment_strlen+1) );
00566
00567 for (i = 0; i < alignment_strlen;i++) {
00568 template_str[i] = '-';
00569 target_str[i] = '-';
00570 template_ss_str[i] = ' ';
00571 target_ss_str[i] = ' ';
00572 align_str[i] = ' ';
00573 core_str[i] = ' ';
00574 }
00575
00576
00577 k = -1;
00578 p = 0;
00579 for (i = 0; i < target_data->len; i++) {
00580 if ( this->target_align[i] == -1) {
00581 target_str[p] = target_data->residue[i].RS;
00582 target_ss_str[p] = target_data->residue[i].SS;
00583 p++;
00584 } else {
00585 for (j = k+1; j < this->target_align[i]; j++) {
00586 template_str[p] = template_data->residue[ j ].RS;
00587 template_ss_str[p] = template_data->residue[ j ].SS;
00588 if (template_data->residue[j].core_num != -1) {
00589 core_str[p] = '-';
00590 }
00591 p++;
00592 }
00593 k = this->target_align[i];
00594 target_str[p] = target_data->residue[i].RS;
00595 target_ss_str[p] = target_data->residue[i].SS;
00596 template_str[p] = template_data->residue[ k ].RS;
00597 template_ss_str[p] = template_data->residue[ k ].SS;
00598 if (template_data->residue[j].core_num != -1) {
00599 core_str[p] = '-';
00600 }
00601 int aa1 = AA1toAANum( template_data->residue[k].RS );
00602 int aa2 = AA1toAANum( target_data->residue[i].RS );
00603 if ( aa1 >= 0 && aa1 < 20 && aa2 >= 0 && aa2 < 20 ) {
00604 align_str[p] = param->ALIGN_MATRIX[ aa1 ][ aa2 ];
00605 }
00606
00607
00608 p++;
00609 }
00610 assert(p < alignment_strlen);
00611 }
00612 for (j = k+1; j < template_data->len; j++) {
00613 template_str[p] = template_data->residue[ j ].RS;
00614 template_ss_str[p] = template_data->residue[ j ].SS;
00615 if (template_data->residue[j].core_num != -1) {
00616 core_str[p] = '-';
00617 }
00618 p++;
00619 }
00620
00621 template_str[p] = 0;
00622 target_str[p] = 0;
00623 target_ss_str[p] = 0;
00624 template_ss_str[p] = 0;
00625 align_str[p] = 0;
00626 core_str[p] = 0;
00627
00628 for (i = 0 ; i < alignment_strlen; i++) {
00629 if (target_ss_str[i] == 'C' || target_ss_str[i] == 'L' )
00630 target_ss_str[i] = ' ';
00631 if (template_ss_str[i] == 'C' || template_ss_str[i] == 'L' )
00632 template_ss_str[i] = ' ';
00633
00634 }
00635
00636 #else
00637 GetAlignmentStrings( param,
00638 template_data,
00639 target_data,
00640 &template_str, &template_ss_str,
00641 &target_str, &target_ss_str,
00642 &align_str, &core_str);
00643 #endif
00644
00645
00646
00647
00648
00649
00650
00651
00652
00653
00654
00655
00656
00657
00658
00659
00660
00661
00662
00663
00664
00665
00666 AlignFeatures features;
00667 features.CalcFeatures( template_data, target_data, this );
00668 sprintf(buffer, "<features>\n");
00669 append_to_str(&outxml, buffer );
00670 for (int i = 0; i < features.Count(); i++) {
00671 if ( features.IsSet(i) ) {
00672 sprintf(buffer, "\t<%s>%g</%s>\n", features.GetName(i), features.Get(i), features.GetName(i));
00673 append_to_str(&outxml, buffer );
00674 }
00675 }
00676 sprintf(buffer, "</features>\n");
00677 append_to_str(&outxml, buffer );
00678
00679
00680
00681 sprintf(buffer, "<alignment >\n<query_ss_>");
00682 append_to_str(&outxml, buffer );
00683 append_to_str(&outxml, target_ss_str);
00684 append_to_str(&outxml, "</query_ss_>\n<query_seq>");
00685 append_to_str(&outxml, target_str);
00686 append_to_str(&outxml, "</query_seq>\n<match_str>");
00687 append_to_str(&outxml, align_str);
00688 append_to_str(&outxml, "</match_str>\n<templ_seq>");
00689 append_to_str(&outxml, template_str );
00690 append_to_str(&outxml, "</templ_seq>\n<templ_ss_>");
00691 append_to_str(&outxml, template_ss_str);
00692 append_to_str(&outxml, "</templ_ss_>\n<core_str_>");
00693 append_to_str(&outxml, core_str);
00694 append_to_str(&outxml, "</core_str_>\n</alignment>\n");
00695
00696 free(target_str);
00697 free(template_str);
00698 free(target_ss_str);
00699 free(align_str);
00700 free(template_ss_str);
00701
00702
00703 if (param->output_template_loc) {
00704 sprintf(buffer, "<template_loc><start><resSeq>%d</resSeq><iCode>%c</iCode></start>\n", template_data->start_res, template_data->start_icode);
00705 append_to_str(&outxml, buffer);
00706 sprintf(buffer, "<stop><resSeq>%d</resSeq><iCode>%c</iCode></stop></template_loc>\n", template_data->stop_res, template_data->stop_icode);
00707 append_to_str(&outxml, buffer);
00708 }
00709
00710
00711 if (param->output_backbone && param->rotamer_array != NULL) {
00712 if (param->v_level > 0)
00713 fprintf(stderr, "Doing Sidechain packing\n");
00714 target_data->CopyBackBone( template_data, this);
00715 target_data->SideChain_Init( param );
00716 target_data->SideChain_Pack( param );
00717 char *pdb_text = target_data->WritePDBText();
00718 sprintf(buffer, "<sidechain_pack><lret>%g</lret><atom_count>%d</atom_count><def_res>%d</def_res><undef_res>%d</undef_res></sidechain_pack>\n",
00719 target_data->Calc_LRET( param ),
00720 target_data->Count_Atoms( ),
00721 target_data->Count_SetResidues( ),
00722 target_data->Count_UnSetResidues( )
00723 );
00724 append_to_str(&outxml, buffer);
00725 append_to_str(&outxml, "<pdb sidechain='yes'>\n");
00726 append_to_str(&outxml, pdb_text);
00727 append_to_str(&outxml, "</pdb>\n");
00728 } else if ( param->output_pdb ) {
00729 target_data->CopyBackBone( template_data, this);
00730 char *pdb_text = target_data->WritePDBText();
00731 append_to_str(&outxml, "<pdb sidechain='no'>\n");
00732 append_to_str(&outxml, pdb_text);
00733 append_to_str(&outxml, "</pdb>\n");
00734 free(pdb_text);
00735 }
00736
00737 sprintf(buffer, "</threading>");
00738 append_to_str(&outxml, buffer);
00739
00740
00741 return outxml;
00742 }
00743
00744
00745
00746 void AlignmentStruct::GetAlignmentStrings( ProspectParam *param,
00747 TemplateStruct *template_data,
00748 TargetStruct *target_data,
00749 char** template_str, char** template_ss_str,
00750 char** target_str, char** target_ss_str,
00751 char** align_str, char** core_str) {
00752
00753 int alignment_strlen = target_data->len + template_data->len;
00754
00755 if (template_str != NULL)
00756 *template_str = (char *)malloc(sizeof(char) * (alignment_strlen+1) );
00757 if (template_ss_str != NULL)
00758 *template_ss_str = (char *)malloc(sizeof(char) * (alignment_strlen+1) );
00759 if (target_str != NULL)
00760 *target_str = (char *)malloc(sizeof(char) * (alignment_strlen+1) );
00761 if (target_ss_str != NULL)
00762 *target_ss_str = (char *)malloc(sizeof(char) * (alignment_strlen+1) );
00763 if ( align_str!= NULL)
00764 *align_str = (char *)malloc(sizeof(char) * (alignment_strlen+1) );
00765 if ( core_str!= NULL)
00766 *core_str = (char *)malloc(sizeof(char) * (alignment_strlen+1) );
00767
00768 for (int i = 0; i < alignment_strlen;i++) {
00769 if ( template_str != NULL)
00770 (*template_str)[i] = '-';
00771 if ( target_str != NULL)
00772 (*target_str)[i] = '-';
00773 if ( template_ss_str != NULL)
00774 (*template_ss_str)[i] = ' ';
00775 if ( target_ss_str != NULL)
00776 (*target_ss_str)[i] = ' ';
00777 if ( align_str != NULL)
00778 (*align_str)[i] = ' ';
00779 if ( core_str != NULL)
00780 (*core_str)[i] = ' ';
00781 }
00782
00783
00784 int k = -1;
00785 int p = 0;
00786 for (int i = 0; i < target_data->len; i++) {
00787 if ( this->target_align[i] == -1) {
00788 if ( target_str != NULL )
00789 (*target_str)[p] = target_data->residue[i].RS;
00790 if ( target_ss_str != NULL)
00791 (*target_ss_str)[p] = target_data->residue[i].SS;
00792 p++;
00793 } else {
00794 int j;
00795 for (j = k+1; j < this->target_align[i]; j++) {
00796 if ( template_str != NULL)
00797 (*template_str)[p] = template_data->residue[ j ].RS;
00798 if ( template_ss_str != NULL )
00799 (*template_ss_str)[p] = template_data->residue[ j ].SS;
00800 if ( core_str != NULL ) {
00801 if (template_data->t_residue_info[j].core_num != -1) {
00802 (*core_str)[p] = '-';
00803 }
00804 }
00805 p++;
00806 }
00807 k = this->target_align[i];
00808 if ( target_str != NULL )
00809 (*target_str)[p] = target_data->residue[i].RS;
00810 if (target_ss_str != NULL)
00811 (*target_ss_str)[p] = target_data->residue[i].SS;
00812 if (template_str != NULL)
00813 (*template_str)[p] = template_data->residue[ k ].RS;
00814 if (template_ss_str != NULL)
00815 (*template_ss_str)[p] = template_data->residue[ k ].SS;
00816 if (core_str != NULL) {
00817 if (template_data->t_residue_info[j].core_num != -1) {
00818 (*core_str)[p] = '-';
00819 }
00820 }
00821 if (align_str != NULL) {
00822 int aa1 = AA1toAANum( template_data->residue[k].RS );
00823 int aa2 = AA1toAANum( target_data->residue[i].RS );
00824 if ( aa1 >= 0 && aa1 < 20 && aa2 >= 0 && aa2 < 20 ) {
00825 if (param != NULL)
00826 (*align_str)[p] = param->ALIGN_MATRIX[ aa1 ][ aa2 ];
00827 else if (aa1 == aa2)
00828 (*align_str)[p] = '|';
00829 }
00830 }
00831 p++;
00832 }
00833 assert(p < alignment_strlen);
00834 }
00835 for (int j = k+1; j < template_data->len; j++) {
00836 if ( template_str != NULL)
00837 (*template_str)[p] = template_data->residue[ j ].RS;
00838 if (template_ss_str != NULL)
00839 (*template_ss_str)[p] = template_data->residue[ j ].SS;
00840 if (core_str != NULL) {
00841 if (template_data->t_residue_info[j].core_num != -1) {
00842 (*core_str)[p] = '-';
00843 }
00844 }
00845 p++;
00846 }
00847 if ( template_str != NULL)
00848 (*template_str)[p] = 0;
00849 if (target_str != NULL)
00850 (*target_str)[p] = 0;
00851 if (target_ss_str != NULL)
00852 (*target_ss_str)[p] = 0;
00853 if (template_ss_str != NULL)
00854 (*template_ss_str)[p] = 0;
00855 if (align_str != NULL)
00856 (*align_str)[p] = 0;
00857 if (core_str != NULL)
00858 (*core_str)[p] = 0;
00859
00860 for (int i = 0 ; i < alignment_strlen; i++) {
00861 if (target_ss_str != NULL && ((*target_ss_str)[i] == 'C' || (*target_ss_str)[i] == 'L') )
00862 (*target_ss_str)[i] = ' ';
00863 if (template_ss_str != NULL && ((*template_ss_str)[i] == 'C' || (*template_ss_str)[i] == 'L') )
00864 (*template_ss_str)[i] = ' ';
00865 }
00866 }
00867
00868 char * AlignmentStruct::GetTargetAlignStr(TemplateStruct *template_data, TargetStruct *target_data) {
00869 char *target_align;
00870 GetAlignmentStrings( NULL, template_data, target_data, NULL, NULL, &target_align, NULL, NULL, NULL);
00871 return target_align;
00872 }
00873
00874
00875 char * AlignmentStruct::GetTemplateAlignStr(TemplateStruct *template_data, TargetStruct *target_data) {
00876 char *template_align;
00877 GetAlignmentStrings( NULL, template_data, target_data, &template_align, NULL, NULL, NULL, NULL, NULL);
00878 return template_align;
00879 }
00880
00881
00882 int AlignmentStruct::SetAlign( AlignmentStruct &rhs, int target_start, int target_end) {
00883 Copy( rhs );
00884 if ( target_start < 0 )
00885 target_align = 0;
00886 if (target_end < 0)
00887 target_end = target_len;
00888 for ( int i = 0; i < target_start; i++)
00889 target_align[i] = -1;
00890 for ( int i = target_end; i < target_len; i++)
00891 target_align[i] = -1;
00892 for ( int i =0; i < template_len; i++)
00893 template_align[i] = -1;
00894 for ( int i = target_start; i < target_end; i++) {
00895 if ( target_align[i] != -1 )
00896 template_align[ target_align[i] ] = i;
00897 }
00898 return 0;
00899 }
00900
00901 int AlignmentStruct::SetAlign( char *target_align_str, char *template_align_str ) {
00902 int target_align_len = 0,
00903 template_align_len = 0;
00904
00905 for (target_align_len = 0;
00906 target_align_str[target_align_len] != '\r' &&
00907 target_align_str[target_align_len] != 0 &&
00908 target_align_str[target_align_len] != '\n'; target_align_len++) { }
00909
00910 for (template_align_len = 0;
00911 template_align_str[template_align_len] != '\r' &&
00912 template_align_str[template_align_len] != 0 &&
00913 template_align_str[template_align_len] != '\n'; template_align_len++) { }
00914
00915 if ( target_align_len != template_align_len )
00916 return 1;
00917 int count_1 = 0;
00918 int count_2 = 0;
00919 target_align = (int *)malloc(sizeof(int) * target_align_len);
00920 template_align = (int *)malloc(sizeof(int) * target_align_len);
00921
00922 for (int i = 0; i < target_align_len; i++) {
00923 target_align[i] = -1;
00924 template_align[i] = -1;
00925 }
00926
00927 for (int i = 0; i < target_align_len; i++) {
00928 if ( target_align_str[i] != '-' && template_align_str[i] != '-' ) {
00929 target_align[count_1] = count_2;
00930 template_align[count_2] = count_1;
00931 }
00932 if ( target_align_str[i] != '-') {
00933 count_1++;
00934 }
00935 if ( template_align_str[i] != '-') {
00936 count_2++;
00937 }
00938 }
00939 target_len = count_1;
00940 template_len = count_2;
00941
00942
00943
00944
00945
00946
00947
00948
00949
00950
00951 return 0;
00952 }
00953
00954
00955 int AlignmentStruct::SetAlignFromNumStr( char *target_align_str, char *template_align_str ) {
00956 target_len = 0;
00957 template_len = 0;
00958 for ( int i = 0; i < strlen( target_align_str ); i++ ) {
00959 if ( target_align_str[i] == ',' )
00960 target_len++;
00961 }
00962 for ( int i = 0; i < strlen( template_align_str ); i++ ) {
00963 if ( template_align_str[i] == ',' )
00964 template_len++;
00965 }
00966
00967 target_align = (int *)malloc(sizeof(int) * target_len );
00968 template_align = (int *)malloc(sizeof(int) * template_len );
00969 target_len = 0;
00970 template_len = 0;
00971
00972 char *ptr;
00973 char *tmp_target = strdup( target_align_str );
00974 ptr = strtok( tmp_target, " ,\t\n" );
00975 if (!ptr)
00976 return 1;
00977 do {
00978 target_align[ target_len ] = atoi( ptr );
00979 target_len++;
00980 } while ( (ptr = strtok( NULL, " ,\t\n" )) != NULL );
00981 free( tmp_target );
00982
00983 char *tmp_template = strdup( template_align_str );
00984 ptr = strtok( tmp_template, " ,\t\n" );
00985 if (!ptr)
00986 return 1;
00987 do {
00988 template_align[ template_len ] = atoi( ptr );
00989 template_len++;
00990 } while ( (ptr = strtok( NULL, " ,\t\n" )) != NULL );
00991 free( tmp_template );
00992
00993 return 0;
00994 }
00995
00996
00997
00998
00999 typedef struct {
01000 int target_start, target_end;
01001 int template_start, template_end;
01002 } AlignPairPos;
01003
01004 int AlignmentStruct::SetAlign( ProspectParam *param,
01005 WeightArray weights,
01006 TemplateStruct *template_data,
01007 TargetStruct *target_data,
01008 int *core_pos, char *core_del ) {
01009
01010
01011 target_align = (int *)malloc(sizeof(int) * target_data->len);
01012 for (int i = 0; i < target_data->len; i++)
01013 target_align[i] = -1;
01014 template_align = (int *)malloc(sizeof(int) * template_data->len);
01015 for (int i = 0; i < template_data->len; i++)
01016 template_align[i] = -1;
01017
01018 target_len = target_data->len;
01019 template_len = template_data->len;
01020
01021
01022 #if 1
01023
01024 ProspectDynamicAlign dynamic_align;
01025 dynamic_align.CalcEnergyTable(param,
01026 weights,
01027 template_data,
01028 target_data);
01029 dynamic_align.Align( );
01030 for (int i = 0; i < template_data->core_count; i++) {
01031 if ( core_pos[i] >= 0 && core_pos[i] < target_data->len ) {
01032 if ( core_del && core_del[i] )
01033 dynamic_align.ForceYDelete( template_data->cores[i].head, template_data->cores[i].len );
01034 else
01035 dynamic_align.ForceAlign( core_pos[i], template_data->cores[i].head, template_data->cores[i].len );
01036 }
01037
01038 if ( core_pos[i] < 0 ) {
01039 dynamic_align.ForceYDelete( 0, template_data->cores[i].head + template_data->cores[i].len );
01040 }
01041 if ( core_pos[i] > target_data->len ) {
01042 dynamic_align.ForceYDelete( template_data->cores[i].head, template_data->len - template_data->cores[i].head );
01043 }
01044
01045 int target_start = core_pos[i] + template_data->cores[ i ].len;
01046 if ( core_del && core_del[i] )
01047 target_start = core_pos[i];
01048 int template_start = template_data->cores[ i ].head + template_data->cores[ i ].len;
01049 if ( target_start >= 0 && target_start < target_data->len && template_start < template_data->len)
01050 dynamic_align.ContinueAlign( target_start , template_start );
01051 }
01052
01053 dynamic_align.TraceAlign( target_align, template_align );
01054 if (param->v_level > 0) {
01055 for (int i = 0; i < target_data->len; i++)
01056 printf("%d ,", target_align[i]);
01057 printf("\n");
01058 for (int i = 0; i < template_data->len; i++)
01059 printf("%d ,", template_align[i]);
01060 printf("\n");
01061 }
01062 if ( param->v_level > 0) {
01063 EnergyArray tmp_energy;
01064 GetEnergies( param, template_data, target_data, tmp_energy);
01065 pValType total = 0;
01066 for (long j = 0; j < e_count; j++)
01067 total += (double)weights[j] * (double)tmp_energy[j];
01068 printf("Alignment Total: %g \n", total);
01069 }
01070
01071
01072
01073 #else
01074
01076
01078
01079
01080
01081 for (int i = 0; i < template_data->core_count; i++) {
01082 if ( core_del == NULL || !core_del[i] ) {
01083 if ( param->v_level > 1 ) {
01084 printf("align; %d,%d -- %d,%d\n",
01085 core_pos[i], template_data->cores[i].head ,
01086 core_pos[i]+ template_data->cores[i].len, template_data->cores[i].head + template_data->cores[i].len );
01087 }
01088 int start_x = core_pos[i];
01089 int start_y = template_data->cores[i].head;
01090 for (int j = 0; j < template_data->cores[i].len; j++) {
01091 int x = start_x + j;
01092 int y = start_y + j;
01093 if ( x >= 0 && x < target_data->len && y >= 0 && y < template_data->len) {
01094 target_align[x] = y;
01095 template_align[y] = x;
01096 }
01097 }
01098 }
01099 }
01100
01101 ProspectDynamicAlign dynamic_align;
01102 dynamic_align.CalcEnergyTable(param,
01103 weights,
01104 template_data,
01105 target_data);
01106
01107 int align_count = 0;
01108
01109 AlignPairPos *align_pos = new AlignPairPos[template_data->core_count + 1];
01110
01111
01112 for ( int i = 0; i < template_data->core_count - 1; i++ ) {
01113 align_pos[align_count].target_start = core_pos[i] + template_data->cores[ i ].len;
01114 if ( core_del && core_del[i] )
01115 align_pos[align_count].target_start = core_pos[i];
01116 align_pos[align_count].target_end = core_pos[i+1];
01117 align_pos[align_count].template_start = template_data->cores[ i ].head + template_data->cores[ i ].len;
01118 align_pos[align_count].template_end = template_data->cores[ i+1 ].head;
01119 if ( align_pos[align_count].target_end > 0 && align_pos[align_count].template_start < template_data->len )
01120 align_count++;
01121 }
01122
01123 if ( template_data->core_count > 0 ) {
01124
01125 align_pos[align_count].target_start = 0;
01126 align_pos[align_count].target_end = core_pos[0];
01127 align_pos[align_count].template_start = 0;
01128 align_pos[align_count].template_end = template_data->cores[ 0 ].head;
01129 if ( align_pos[align_count].target_end > 0 && align_pos[align_count].template_start < template_data->len )
01130 align_count++;
01131
01132 align_pos[align_count].target_start = core_pos[ template_data->core_count - 1 ] +
01133 template_data->cores[ template_data->core_count - 1 ].len;
01134 if ( core_del && core_del[ template_data->core_count - 1 ] )
01135 align_pos[align_count].target_start = core_pos[ template_data->core_count - 1 ];
01136 align_pos[align_count].target_end = target_data->len;
01137 align_pos[align_count].template_start = template_data->cores[ template_data->core_count - 1 ].head +
01138 template_data->cores[ template_data->core_count - 1 ].len ;
01139 align_pos[align_count].template_end = template_data->len;
01140 if ( align_pos[align_count].target_end > 0 && align_pos[align_count].template_start < template_data->len )
01141 align_count++;
01142 }
01143
01144 for ( int i = 0; i < align_count; i++ ) {
01145 if ( align_pos[i].target_start < 0 )
01146 align_pos[i].target_start = 0;
01147 if ( align_pos[i].target_end < 0 )
01148 align_pos[i].target_end = 0;
01149 if ( align_pos[i].target_start > target_data->len - 1 )
01150 align_pos[i].target_start = target_data->len - 1;
01151 if ( align_pos[i].target_end > target_data->len )
01152 align_pos[i].target_end = target_data->len;
01153 if ( align_pos[i].template_start > template_data->len - 1 )
01154 align_pos[i].template_start = template_data->len - 1;
01155 if ( align_pos[i].template_end > template_data->len )
01156 align_pos[i].template_end = template_data->len;
01157 }
01158
01159 for ( int i = 0; i < align_count; i++ ) {
01160 if (param->v_level > 1 ) {
01161 printf("align: %d,%d -- %d,%d\n",
01162 align_pos[i].target_start, align_pos[i].template_start,
01163 align_pos[i].target_end, align_pos[i].template_end );
01164 }
01165 dynamic_align.Align( align_pos[i].target_start, align_pos[i].template_start,
01166 align_pos[i].target_end, align_pos[i].template_end );
01167 dynamic_align.TraceAlign(align_pos[i].target_start, align_pos[i].template_start,
01168 align_pos[i].target_end, align_pos[i].template_end,
01169 target_align, template_align );
01170 }
01171 delete align_pos;
01172 if (param->v_level > 0) {
01173 for (int i = 0; i < target_data->len; i++)
01174 printf("%d ,", target_align[i]);
01175 printf("\n");
01176 for (int i = 0; i < template_data->len; i++)
01177 printf("%d ,", template_align[i]);
01178 printf("\n");
01179 }
01180 if ( param->v_level > 0) {
01181 EnergyArray tmp_energy;
01182 GetEnergies( param, template_data, target_data, tmp_energy);
01183 pValType total = 0;
01184 for (long j = 0; j < e_count; j++)
01185 total += (double)weights[j] * (double)tmp_energy[j];
01186 printf("Alignment Total: %f \n", total);
01187 }
01188 #endif
01189 return 0;
01190 }
