00001
00002 #include <stdio.h>
00003 #include <stdlib.h>
00004 #include <string.h>
00005 #include <math.h>
00006 #include <assert.h>
00007 #include <string>
00008
00009 #include "open_prospect.h"
00010 #include "tree_decomp.h"
00011 #include "sm_matrix.h"
00012
00013
00014
00015
00016
00017
00018 TargetAlignmentStruct::TargetAlignmentStruct() {
00019 target_1_align = NULL;
00020 target_2_align = NULL;
00021 target_1_len = 0;
00022 target_2_len = 0;
00023 }
00024
00025
00026 TargetAlignmentStruct::TargetAlignmentStruct(const TargetAlignmentStruct &rhs) {
00027 target_1_align = NULL;
00028 target_2_align = NULL;
00029 target_1_len = 0;
00030 target_2_len = 0;
00031 if (rhs.target_1_len > 0) {
00032 target_1_len = rhs.target_1_len;
00033 target_1_align = (int *)malloc(sizeof(int) * target_1_len);
00034 memcpy( target_1_align, rhs.target_1_align, sizeof(int) * target_1_len );
00035 }
00036 if (rhs.target_2_len > 0) {
00037 target_2_len = rhs.target_1_len;
00038 target_2_align = (int *)malloc(sizeof(int) * target_2_len);
00039 memcpy( target_2_align, rhs.target_2_align, sizeof(int) * target_2_len);
00040 }
00041 }
00042
00043 TargetAlignmentStruct::TargetAlignmentStruct(const AlignmentStruct &rhs) {
00044 target_1_align = NULL;
00045 target_2_align = NULL;
00046 target_1_len = 0;
00047 target_2_len = 0;
00048 if (rhs.target_len > 0) {
00049 target_1_len = rhs.target_len;
00050 target_1_align = (int *)malloc(sizeof(int) * target_1_len);
00051 memcpy( target_1_align, rhs.target_align, sizeof(int) * target_1_len );
00052 }
00053 if (rhs.template_len > 0) {
00054 target_2_len = rhs.template_len;
00055 target_2_align = (int *)malloc(sizeof(int) * target_2_len);
00056 memcpy( target_2_align, rhs.template_align, sizeof(int) * target_2_len);
00057 }
00058 }
00059
00060
00061
00062 TargetAlignmentStruct::~TargetAlignmentStruct() {
00063 if (target_1_align != NULL)
00064 free(target_1_align);
00065 if (target_2_align != NULL)
00066 free(target_2_align);
00067 }
00068
00069
00070 void TargetAlignmentStruct::SetAlign(TargetStruct *target_1_in, int *target_1_align_in,
00071 TargetStruct *target_2_in, int *target_2_align_in) {
00072 target_1_len = target_1_in->len;
00073 target_2_len = target_2_in->len;
00074 target_1_align = target_1_align_in;
00075 target_2_align = target_2_align_in;
00076 }
00077
00078
00079 void TargetAlignmentStruct::SetAlign(TargetStruct *target_1_in, int target_1_offset,
00080 TargetStruct *target_2_in, int target_2_offset,
00081 int len) {
00082 target_1_len = target_1_in->len;
00083 target_2_len = target_2_in->len;
00084 target_1_align = (int *)malloc(sizeof(int) * target_1_len);
00085 target_2_align = (int *)malloc(sizeof(int) * target_2_len);
00086 ClearAlign();
00087 for (int i = 0; i < len; i++) {
00088 if ( target_1_offset + i < target_1_len && target_2_offset + i < target_2_len ) {
00089 target_1_align[ target_1_offset + i ] = target_2_offset + i;
00090 target_2_align[ target_2_offset + i ] = target_1_offset + i;
00091 }
00092 }
00093 }
00094
00095
00096 void TargetAlignmentStruct::SetAlign(int *target_1, int *target_2, int len_both) {
00097 target_1_align = target_1;
00098 target_2_align = target_2;
00099 target_1_len = len_both;
00100 target_2_len = len_both;
00101 }
00102
00103
00104 void TargetAlignmentStruct::SetAlign(int *target_1, int len_1, int *target_2, int len_2) {
00105 target_1_align = target_1;
00106 target_2_align = target_2;
00107 target_1_len = len_1;
00108 target_2_len = len_2;
00109 }
00110
00111
00112 char TargetAlignmentStruct::SetAlign(const char *target1, const char *target2) {
00113 target_1_len = 0;
00114 for ( int i = 0; target1[i] != 0 && target1[i] != '\n'; i++) {
00115 int aa_num = AA1toAANum( target1[i] );
00116 if ( aa_num >= 0 && aa_num < 20 )
00117 target_1_len++;
00118 }
00119 target_2_len = 0;
00120 for ( int i = 0; target2[i] != 0 && target2[i] != '\n'; i++) {
00121 int aa_num = AA1toAANum( target2[i] );
00122 if ( aa_num >= 0 && aa_num < 20 )
00123 target_2_len++;
00124 }
00125 if ( target_1_len == 0 || target_2_len == 0 )
00126 return 1;
00127
00128 target_1_align = (int *)malloc(sizeof(int) * target_1_len);
00129 target_2_align = (int *)malloc(sizeof(int) * target_2_len);
00130 for (int i = 0; i < target_1_len; i++)
00131 target_1_align[i] = -1;
00132 for (int i = 0; i < target_2_len; i++)
00133 target_2_align[i] = -1;
00134 int count_1 = 0, count_2 = 0;
00135 for (int i = 0; target1[i] != 0 && target1[i] != '\n'; i++) {
00136 int aa_num1 = AA1toAANum( target1[i] );
00137 int aa_num2 = AA1toAANum( target2[i] );
00138 if ( ( aa_num1 >= 0 && aa_num1 < 20 ) && ( aa_num2 >= 0 && aa_num2 < 20 ) ) {
00139 target_1_align[count_1] = count_2;
00140 target_2_align[count_2] = count_1;
00141 }
00142 if ( aa_num1 >= 0 && aa_num1 < 20 ) {
00143 count_1++;
00144 }
00145 if ( aa_num2 >= 0 && aa_num2 < 20 ) {
00146 count_2++;
00147 }
00148 }
00149 return CheckAlign();
00150 }
00151
00152 char TargetAlignmentStruct::LoadFast(char *fast_file) {
00153 std::string string1, string2;
00154 char buffer[3000];
00155 FILE *file = fopen( fast_file, "r" );
00156 if (file == NULL)
00157 return 1;
00158 while (fgets( buffer, sizeof(buffer), file) ) {
00159 char *ptr;
00160 if ( NULL != (ptr = strstr( buffer, "1:" )) ) {
00161 ptr += 3;
00162 for (int i = 0, j = 0; ptr[i] != 0; i++) {
00163 if ( ptr[i] != '\n' && ptr[i] != ' ' && ptr[i] != '*' ) {
00164 ptr[j] = ptr[i];
00165 j++;
00166 } else
00167 ptr[j] = 0;
00168 }
00169 string1 += ptr;
00170 } else if ( NULL != (ptr = strstr( buffer, "2:" )) ) {
00171 ptr += 3;
00172 for (int i = 0, j = 0; ptr[i] != 0; i++) {
00173 if ( ptr[i] != '\n' && ptr[i] != ' ' && ptr[i] != '*' ) {
00174 ptr[j] = ptr[i];
00175 j++;
00176 } else
00177 ptr[j] = 0;
00178 }
00179 string2 += ptr;
00180 }
00181 }
00182
00183 fclose(file);
00184
00185 return SetAlign( string1.c_str(), string2.c_str() );
00186 }
00187
00188
00189 char TargetAlignmentStruct::MergeAlign( TargetAlignmentStruct *align_1, TargetAlignmentStruct *align_2) {
00190
00191 if (align_1->target_1_len != align_2->target_1_len || align_1->target_2_len != align_2->target_2_len )
00192 return 1;
00193
00194 target_1_len = align_1->target_1_len;
00195 target_2_len = align_1->target_2_len;
00196 target_1_align = (int *)malloc(sizeof(int) * target_1_len);
00197 target_2_align = (int *)malloc(sizeof(int) * target_2_len);
00198 ClearAlign();
00199 for (int i = 0; i < target_1_len; i++) {
00200 if (align_1->target_1_align[i] != -1) {
00201 target_1_align[i] = align_1->target_1_align[i];
00202 target_2_align[ align_1->target_1_align[i] ] = i;
00203 }
00204 if (align_2->target_1_align[i] != -1) {
00205 target_1_align[i] = align_2->target_1_align[i];
00206 target_2_align[ align_2->target_1_align[i] ] = i;
00207 }
00208 }
00209 return CheckAlign();
00210 }
00211
00212 void TargetAlignmentStruct::ClearAlign(void) {
00213 for ( int i = 0; i < target_1_len; i++) {
00214 target_1_align[i] = -1;
00215 }
00216 for ( int i = 0; i < target_2_len; i++) {
00217 target_2_align[i] = -1;
00218 }
00219 }
00220
00221 void TargetAlignmentStruct::ModAlign(int target_1_offset, int target_2_offset, int len) {
00222 for (int i = 0; i < len; i++) {
00223 if ( (target_1_offset + i) >= 0 && (target_1_offset + i) < target_1_len &&
00224 (target_2_offset + i) >= 0 && (target_1_offset + i) < target_2_len) {
00225 target_1_align[ target_1_offset + i ] = target_2_offset + i;
00226 target_2_align[ target_2_offset + i ] = target_1_offset + i;
00227 }
00228 }
00229 }
00230
00231 char TargetAlignmentStruct::CheckAlign(void) {
00232 int last_pos = -1;
00233 for ( int i = 0; i < target_1_len; i++) {
00234 if ( target_1_align[i] > last_pos ) {
00235 last_pos = target_1_align[i];
00236 } else if (target_1_align[i] != -1) {
00237 return 1;
00238 }
00239 }
00240 last_pos = -1;
00241 for ( int i = 0; i < target_2_len; i++) {
00242 if ( target_2_align[i] > last_pos ) {
00243 last_pos = target_2_align[i];
00244 } else if (target_2_align[i] != -1) {
00245 return 1;
00246 }
00247 }
00248 return 0;
00249 }
00250
00251
00252 pValType TargetAlignmentStruct::GetCaRmsd(TargetStruct *target_1, TargetStruct *target_2) {
00253 vec3 trans1, trans2;
00254 mat3 rot;
00255 return Get1to2Trans( target_1, target_2, trans1, rot, trans2 );
00256 }
00257
00258
00259 pValType TargetAlignmentStruct::AlignTargetCa(TargetStruct *target_1, TargetStruct *target_2,
00260 int *target_1_align, int *target_2_align,
00261 vec3 trans1, mat3 rot, vec3 trans2) {
00262
00263 int align_count = 0;
00264 for (int i = 0; i < target_1->len; i++) {
00265 if ( target_1_align[i] != -1 )
00266 align_count++;
00267 }
00268
00269 vec3 ca_vec_1[ align_count ], ca_vec_2[ align_count ];
00270 align_count = 0;
00271 for (int i = 0; i < target_1->len; i++) {
00272 if ( target_1_align[i] != -1 ) {
00273 ca_vec_1[align_count][0] = target_1->structure[ i ].atom[ (int)target_1->structure[i].bb_num[ anCa ] ][0];
00274 ca_vec_1[align_count][1] = target_1->structure[ i ].atom[ (int)target_1->structure[i].bb_num[ anCa ] ][1];
00275 ca_vec_1[align_count][2] = target_1->structure[ i ].atom[ (int)target_1->structure[i].bb_num[ anCa ] ][2];
00276 int i_2 = target_1_align[i];
00277 ca_vec_2[align_count][0] = target_2->structure[ i_2 ].atom[ (int)target_2->structure[i_2].bb_num[ anCa ] ][0];
00278 ca_vec_2[align_count][1] = target_2->structure[ i_2 ].atom[ (int)target_2->structure[i_2].bb_num[ anCa ] ][1];
00279 ca_vec_2[align_count][2] = target_2->structure[ i_2 ].atom[ (int)target_2->structure[i_2].bb_num[ anCa ] ][2];
00280 align_count++;
00281 }
00282 }
00283
00284 smCalcMeanOrigin( ca_vec_1, align_count, trans1 );
00285 smCalcMeanOrigin( ca_vec_2, align_count, trans2 );
00286 smMinus( trans2 );
00287
00288 smMoveMeanOrigin( ca_vec_1, align_count );
00289 smMoveMeanOrigin( ca_vec_2, align_count );
00290
00291 if ( !smCalcMinRmsdRotation( ca_vec_1, ca_vec_2, align_count, rot) ) {
00292 smMultMat(rot, ca_vec_1, align_count);
00293 float rmsd = smCalcRmsd( ca_vec_1, ca_vec_2, align_count );
00294 return rmsd;
00295 }
00296 return 1000;
00297 }
00298
00299
00300 pValType TargetAlignmentStruct::AlignVec3( vec3 *target_1, vec3 *target_2, vec3 trans1, mat3 rot, vec3 trans2) {
00301 int align_count = 0;
00302 for (int i = 0; i < target_1_len; i++) {
00303 if ( target_1_align[i] != -1 )
00304 align_count++;
00305 }
00306 vec3 ca_vec_1[ align_count ], ca_vec_2[ align_count ];
00307 align_count = 0;
00308 for (int i = 0; i < target_1_len; i++) {
00309 if ( target_1_align[i] != -1 ) {
00310 ca_vec_1[align_count][0] = target_1[ i ][0];
00311 ca_vec_1[align_count][1] = target_1[ i ][1];
00312 ca_vec_1[align_count][2] = target_1[ i ][2];
00313 int i_2 = target_1_align[i];
00314 ca_vec_2[align_count][0] = target_2[ i_2 ][0];
00315 ca_vec_2[align_count][1] = target_2[ i_2 ][1];
00316 ca_vec_2[align_count][2] = target_2[ i_2 ][2];
00317 align_count++;
00318 }
00319 }
00320 smCalcMeanOrigin( ca_vec_1, align_count, trans1 );
00321 smCalcMeanOrigin( ca_vec_2, align_count, trans2 );
00322 smMinus( trans2 );
00323 smMoveMeanOrigin( ca_vec_1, align_count );
00324 smMoveMeanOrigin( ca_vec_2, align_count );
00325 if ( !smMoveMinRmsd( ca_vec_1, ca_vec_2, align_count ) ) {
00326 float rmsd = smCalcRmsd( ca_vec_1, ca_vec_2, align_count );
00327 return rmsd;
00328 }
00329 return 1000;
00330 }
00331
00332
00333
00334 pValType TargetAlignmentStruct::Get1to2Trans( TargetStruct *target_1, TargetStruct *target_2,
00335 vec3 trans1, mat3 rot, vec3 trans2 ) {
00336 return AlignTargetCa(target_1, target_2, target_1_align, target_2_align,
00337 trans1, rot, trans2);
00338 }
00339
00340 pValType TargetAlignmentStruct::Get2to1Trans( TargetStruct *target_1, TargetStruct *target_2,
00341 vec3 trans1, mat3 rot, vec3 trans2 ) {
00342 return AlignTargetCa(target_2, target_1, target_2_align, target_1_align,
00343 trans1, rot, trans2);
00344 }
00345
00346
00347
00348 void TargetAlignmentStruct::Move1to2(TargetStruct *target_1, TargetStruct *target_2) {
00349 vec3 trans1, trans2;
00350 mat3 rot;
00351
00352 AlignTargetCa(target_1, target_2, target_1_align, target_2_align,
00353 trans1, rot, trans2);
00354 for (int i = 0; i < target_1->len; i++) {
00355 smAdd( trans1, target_1->structure[i].atom, target_1->structure[i].atom_count );
00356 smMultMat( rot, target_1->structure[i].atom, target_1->structure[i].atom_count );
00357 smAdd( trans2, target_1->structure[i].atom, target_1->structure[i].atom_count );
00358 }
00359 }
00360
00361
00362
00363 void TargetAlignmentStruct::Move2to1(TargetStruct *target_1, TargetStruct *target_2) {
00364 vec3 trans1, trans2;
00365 mat3 rot;
00366
00367 AlignTargetCa(target_2, target_1, target_2_align, target_1_align,
00368 trans1, rot, trans2);
00369 for (int i = 0; i < target_2->len; i++) {
00370 smAdd( trans1, target_2->structure[i].atom, target_2->structure[i].atom_count );
00371 smMultMat( rot, target_2->structure[i].atom, target_2->structure[i].atom_count );
00372 smAdd( trans2, target_2->structure[i].atom, target_2->structure[i].atom_count );
00373 }
00374 }
00375
00376
00377 void TargetAlignmentStruct::Move1to2( vec3 *target_1, vec3 *target_2 ) {
00378 vec3 trans1, trans2;
00379 mat3 rot;
00380 AlignVec3(target_1, target_2, trans1, rot, trans2);
00381 smAdd( trans1, target_1, target_1_len );
00382 smMultMat( rot, target_1, target_1_len );
00383 smAdd( trans2, target_1, target_1_len );
00384 }
00385
00386
00387 void TargetAlignmentStruct::ApplyAlignTrans( TargetStruct *target, vec3 trans1, mat3 rot, vec3 trans2) {
00388 ApplyAlignTrans( target->structure, target->len, trans1, rot, trans2 );
00389 }
00390
00391 void TargetAlignmentStruct::ApplyAlignTrans( ResidueAtoms *target, int target_len, vec3 trans1, mat3 rot, vec3 trans2) {
00392 for (int i = 0; i < target_len; i++) {
00393 smAdd( trans1, target[i].atom, target[i].atom_count );
00394 smMultMat( rot, target[i].atom, target[i].atom_count );
00395 smAdd( trans2, target[i].atom, target[i].atom_count );
00396 }
00397 }
00398
00399
00400 void TargetAlignmentStruct::ApplyAlignTrans( vec3 *target, int target_len, vec3 trans1, mat3 rot, vec3 trans2) {
00401 smAdd( trans1, target, target_len );
00402 smMultMat( rot, target, target_len );
00403 smAdd( trans2, target, target_len );
00404 }
00405
