| 1 | #include <stdio.h> |
|---|
| 2 | #include <stdlib.h> |
|---|
| 3 | #include <string.h> |
|---|
| 4 | #include <math.h> |
|---|
| 5 | #include "clustalw.h" |
|---|
| 6 | #define ENDALN 127 |
|---|
| 7 | |
|---|
| 8 | #define MAX(a,b) ((a)>(b)?(a):(b)) |
|---|
| 9 | #define MIN(a,b) ((a)<(b)?(a):(b)) |
|---|
| 10 | |
|---|
| 11 | /* |
|---|
| 12 | * Prototypes |
|---|
| 13 | */ |
|---|
| 14 | static lint pdiff(sint A,sint B,sint i,sint j,sint go1,sint go2); |
|---|
| 15 | static lint prfscore(sint n, sint m); |
|---|
| 16 | static sint gap_penalty1(sint i, sint j,sint k); |
|---|
| 17 | static sint open_penalty1(sint i, sint j); |
|---|
| 18 | static sint ext_penalty1(sint i, sint j); |
|---|
| 19 | static sint gap_penalty2(sint i, sint j,sint k); |
|---|
| 20 | static sint open_penalty2(sint i, sint j); |
|---|
| 21 | static sint ext_penalty2(sint i, sint j); |
|---|
| 22 | static void padd(sint k); |
|---|
| 23 | static void pdel(sint k); |
|---|
| 24 | static void palign(void); |
|---|
| 25 | static void ptracepath(sint *alen); |
|---|
| 26 | static void add_ggaps(void); |
|---|
| 27 | static char * add_ggaps_mask(char *mask, int len, char *path1, char *path2); |
|---|
| 28 | |
|---|
| 29 | /* |
|---|
| 30 | * Global variables |
|---|
| 31 | */ |
|---|
| 32 | extern double **tmat; |
|---|
| 33 | extern float gap_open, gap_extend; |
|---|
| 34 | extern float transition_weight; |
|---|
| 35 | extern sint gap_pos1, gap_pos2; |
|---|
| 36 | extern sint max_aa; |
|---|
| 37 | extern sint nseqs; |
|---|
| 38 | extern sint *seqlen_array; |
|---|
| 39 | extern sint *seq_weight; |
|---|
| 40 | extern sint debug; |
|---|
| 41 | extern Boolean neg_matrix; |
|---|
| 42 | extern sint mat_avscore; |
|---|
| 43 | extern short blosum30mt[], blosum40mt[], blosum45mt[]; |
|---|
| 44 | extern short blosum62mt2[], blosum80mt[]; |
|---|
| 45 | extern short pam20mt[], pam60mt[]; |
|---|
| 46 | extern short pam120mt[], pam160mt[], pam350mt[]; |
|---|
| 47 | extern short gon40mt[], gon80mt[]; |
|---|
| 48 | extern short gon120mt[], gon160mt[], gon250mt[], gon300mt[]; |
|---|
| 49 | extern short clustalvdnamt[],swgapdnamt[]; |
|---|
| 50 | extern short idmat[]; |
|---|
| 51 | extern short usermat[]; |
|---|
| 52 | extern short userdnamat[]; |
|---|
| 53 | extern short def_dna_xref[],def_aa_xref[],dna_xref[],aa_xref[]; |
|---|
| 54 | extern sint max_aln_length; |
|---|
| 55 | extern Boolean distance_tree; |
|---|
| 56 | extern Boolean dnaflag; |
|---|
| 57 | extern char mtrxname[]; |
|---|
| 58 | extern char dnamtrxname[]; |
|---|
| 59 | extern char **seq_array; |
|---|
| 60 | extern char *amino_acid_codes; |
|---|
| 61 | extern char *gap_penalty_mask1,*gap_penalty_mask2; |
|---|
| 62 | extern char *sec_struct_mask1,*sec_struct_mask2; |
|---|
| 63 | extern sint struct_penalties1, struct_penalties2; |
|---|
| 64 | extern Boolean use_ss1, use_ss2; |
|---|
| 65 | |
|---|
| 66 | static sint print_ptr,last_print; |
|---|
| 67 | static sint *displ; |
|---|
| 68 | |
|---|
| 69 | static char **alignment; |
|---|
| 70 | static sint *aln_len; |
|---|
| 71 | static sint *aln_weight; |
|---|
| 72 | static char *aln_path1, *aln_path2; |
|---|
| 73 | static sint alignment_len; |
|---|
| 74 | static sint **profile1, **profile2; |
|---|
| 75 | static lint *HH, *DD, *RR, *SS; |
|---|
| 76 | static lint *gS; |
|---|
| 77 | static sint matrix[NUMRES][NUMRES]; |
|---|
| 78 | static sint nseqs1, nseqs2; |
|---|
| 79 | static sint prf_length1, prf_length2; |
|---|
| 80 | static sint *gaps; |
|---|
| 81 | static sint gapcoef1,gapcoef2; |
|---|
| 82 | static sint lencoef1,lencoef2; |
|---|
| 83 | static Boolean switch_profiles; |
|---|
| 84 | |
|---|
| 85 | lint prfalign(sint *group, sint *aligned) |
|---|
| 86 | { |
|---|
| 87 | |
|---|
| 88 | static sint i, j, count = 0; |
|---|
| 89 | static sint NumSeq; |
|---|
| 90 | static sint len, len1, len2, is, minlen; |
|---|
| 91 | static sint se1, se2, sb1, sb2; |
|---|
| 92 | static sint maxres; |
|---|
| 93 | static sint int_scale; |
|---|
| 94 | static short *matptr; |
|---|
| 95 | static short *mat_xref; |
|---|
| 96 | static char c; |
|---|
| 97 | static lint score; |
|---|
| 98 | static float scale; |
|---|
| 99 | static double logmin,logdiff; |
|---|
| 100 | static double pcid; |
|---|
| 101 | |
|---|
| 102 | |
|---|
| 103 | alignment = (char **) ckalloc( nseqs * sizeof (char *) ); |
|---|
| 104 | aln_len = (sint *) ckalloc( nseqs * sizeof (sint) ); |
|---|
| 105 | aln_weight = (sint *) ckalloc( nseqs * sizeof (sint) ); |
|---|
| 106 | |
|---|
| 107 | for (i=0;i<nseqs;i++) |
|---|
| 108 | if (aligned[i+1] == 0) group[i+1] = 0; |
|---|
| 109 | |
|---|
| 110 | nseqs1 = nseqs2 = 0; |
|---|
| 111 | for (i=0;i<nseqs;i++) |
|---|
| 112 | { |
|---|
| 113 | if (group[i+1] == 1) nseqs1++; |
|---|
| 114 | else if (group[i+1] == 2) nseqs2++; |
|---|
| 115 | } |
|---|
| 116 | |
|---|
| 117 | if ((nseqs1 == 0) || (nseqs2 == 0)) return(0.0); |
|---|
| 118 | |
|---|
| 119 | if (nseqs2 > nseqs1) |
|---|
| 120 | { |
|---|
| 121 | switch_profiles = TRUE; |
|---|
| 122 | for (i=0;i<nseqs;i++) |
|---|
| 123 | { |
|---|
| 124 | if (group[i+1] == 1) group[i+1] = 2; |
|---|
| 125 | else if (group[i+1] == 2) group[i+1] = 1; |
|---|
| 126 | } |
|---|
| 127 | } |
|---|
| 128 | else |
|---|
| 129 | switch_profiles = FALSE; |
|---|
| 130 | |
|---|
| 131 | int_scale = 100; |
|---|
| 132 | if (dnaflag) |
|---|
| 133 | { |
|---|
| 134 | scale=1.0; |
|---|
| 135 | if (strcmp(dnamtrxname, "iub") == 0) |
|---|
| 136 | { |
|---|
| 137 | matptr = swgapdnamt; |
|---|
| 138 | mat_xref = def_dna_xref; |
|---|
| 139 | } |
|---|
| 140 | else if (strcmp(dnamtrxname, "clustalw") == 0) |
|---|
| 141 | { |
|---|
| 142 | matptr = clustalvdnamt; |
|---|
| 143 | mat_xref = def_dna_xref; |
|---|
| 144 | scale=0.66; |
|---|
| 145 | } |
|---|
| 146 | else |
|---|
| 147 | { |
|---|
| 148 | matptr = userdnamat; |
|---|
| 149 | mat_xref = dna_xref; |
|---|
| 150 | } |
|---|
| 151 | maxres = get_matrix(matptr, mat_xref, matrix, neg_matrix, int_scale); |
|---|
| 152 | if (maxres == 0) return((sint)-1); |
|---|
| 153 | matrix[0][4]=transition_weight*matrix[0][0]; |
|---|
| 154 | matrix[4][0]=transition_weight*matrix[0][0]; |
|---|
| 155 | matrix[2][11]=transition_weight*matrix[0][0]; |
|---|
| 156 | matrix[11][2]=transition_weight*matrix[0][0]; |
|---|
| 157 | matrix[2][12]=transition_weight*matrix[0][0]; |
|---|
| 158 | matrix[12][2]=transition_weight*matrix[0][0]; |
|---|
| 159 | } |
|---|
| 160 | else |
|---|
| 161 | { |
|---|
| 162 | |
|---|
| 163 | /* |
|---|
| 164 | calculate the mean of the sequence pc identities between the two groups |
|---|
| 165 | */ |
|---|
| 166 | count = 0; |
|---|
| 167 | pcid = 0.0; |
|---|
| 168 | for (i=0;i<nseqs;i++) |
|---|
| 169 | { |
|---|
| 170 | if (group[i+1] == 1) |
|---|
| 171 | for (j=0;j<nseqs;j++) |
|---|
| 172 | if (group[j+1] == 2) |
|---|
| 173 | { |
|---|
| 174 | count++; |
|---|
| 175 | if (pcid < tmat[i+1][j+1]) pcid = tmat[i+1][j+1]; |
|---|
| 176 | /* |
|---|
| 177 | pcid += tmat[i+1][j+1]; |
|---|
| 178 | */ |
|---|
| 179 | } |
|---|
| 180 | } |
|---|
| 181 | /* |
|---|
| 182 | pcid = pcid/(float)count; |
|---|
| 183 | */ |
|---|
| 184 | if (debug > 0) fprintf(stdout,"mean tmat %3.1f\n", pcid); |
|---|
| 185 | |
|---|
| 186 | scale = 0.75; |
|---|
| 187 | if (strcmp(mtrxname, "blosum") == 0) |
|---|
| 188 | { |
|---|
| 189 | if (distance_tree == FALSE) matptr = blosum40mt; |
|---|
| 190 | else if (pcid > 80.0) |
|---|
| 191 | { |
|---|
| 192 | scale=0.5; |
|---|
| 193 | matptr = blosum80mt; |
|---|
| 194 | scale *= -3.5 + pcid/20.0; |
|---|
| 195 | } |
|---|
| 196 | else if (pcid > 60.0) |
|---|
| 197 | { |
|---|
| 198 | scale=0.5; |
|---|
| 199 | matptr = blosum62mt2; |
|---|
| 200 | scale *= -2.5 + pcid/20.0; |
|---|
| 201 | } |
|---|
| 202 | else if (pcid > 30.0) |
|---|
| 203 | { |
|---|
| 204 | scale=0.5; |
|---|
| 205 | matptr = blosum45mt; |
|---|
| 206 | scale *= -0.4 + pcid/30.0; |
|---|
| 207 | } |
|---|
| 208 | else if (pcid > 10.0) |
|---|
| 209 | { |
|---|
| 210 | matptr = blosum30mt; |
|---|
| 211 | scale *= pcid/30.0; |
|---|
| 212 | } |
|---|
| 213 | else |
|---|
| 214 | { |
|---|
| 215 | matptr = blosum30mt; |
|---|
| 216 | scale *= 0.30; |
|---|
| 217 | } |
|---|
| 218 | mat_xref = def_aa_xref; |
|---|
| 219 | |
|---|
| 220 | } |
|---|
| 221 | else if (strcmp(mtrxname, "pam") == 0) |
|---|
| 222 | { |
|---|
| 223 | if (distance_tree == FALSE) matptr = pam120mt; |
|---|
| 224 | else if (pcid > 80.0) matptr = pam20mt; |
|---|
| 225 | else if (pcid > 60.0) matptr = pam60mt; |
|---|
| 226 | else if (pcid > 40.0) matptr = pam120mt; |
|---|
| 227 | else matptr = pam350mt; |
|---|
| 228 | mat_xref = def_aa_xref; |
|---|
| 229 | } |
|---|
| 230 | else if (strcmp(mtrxname, "gonnet") == 0) |
|---|
| 231 | { |
|---|
| 232 | if (distance_tree == FALSE) matptr = gon120mt; |
|---|
| 233 | else if (pcid > 65.0) |
|---|
| 234 | { |
|---|
| 235 | matptr = gon40mt; |
|---|
| 236 | scale *= -0.5 + pcid/65.0; |
|---|
| 237 | } |
|---|
| 238 | else if (pcid > 45.0) |
|---|
| 239 | { |
|---|
| 240 | matptr = gon80mt; |
|---|
| 241 | scale *= -1.5 + pcid/20.0; |
|---|
| 242 | } |
|---|
| 243 | else if (pcid > 35.0) |
|---|
| 244 | { |
|---|
| 245 | matptr = gon120mt; |
|---|
| 246 | scale *= -2.5 + pcid/10.0; |
|---|
| 247 | } |
|---|
| 248 | else if (pcid > 25.0) |
|---|
| 249 | { |
|---|
| 250 | matptr = gon160mt; |
|---|
| 251 | scale *= -1.5 + pcid/10.0; |
|---|
| 252 | } |
|---|
| 253 | else if (pcid > 15.0) |
|---|
| 254 | { |
|---|
| 255 | matptr = gon250mt; |
|---|
| 256 | scale *= -0.5 + pcid/10.0; |
|---|
| 257 | } |
|---|
| 258 | else |
|---|
| 259 | { |
|---|
| 260 | matptr = gon300mt; |
|---|
| 261 | scale *= 0.5; |
|---|
| 262 | } |
|---|
| 263 | mat_xref = def_aa_xref; |
|---|
| 264 | int_scale /= 10; |
|---|
| 265 | } |
|---|
| 266 | else if (strcmp(mtrxname, "id") == 0) |
|---|
| 267 | { |
|---|
| 268 | matptr = idmat; |
|---|
| 269 | mat_xref = def_aa_xref; |
|---|
| 270 | } |
|---|
| 271 | else |
|---|
| 272 | { |
|---|
| 273 | matptr = usermat; |
|---|
| 274 | mat_xref = aa_xref; |
|---|
| 275 | } |
|---|
| 276 | |
|---|
| 277 | maxres = get_matrix(matptr, mat_xref, matrix, neg_matrix, int_scale); |
|---|
| 278 | if (maxres == 0) |
|---|
| 279 | { |
|---|
| 280 | fprintf(stdout,"Error: matrix %s not found\n", mtrxname); |
|---|
| 281 | return(-1); |
|---|
| 282 | } |
|---|
| 283 | } |
|---|
| 284 | |
|---|
| 285 | /* |
|---|
| 286 | Make the first profile. |
|---|
| 287 | */ |
|---|
| 288 | prf_length1 = 0; |
|---|
| 289 | nseqs1 = 0; |
|---|
| 290 | if (debug>0) fprintf(stdout,"sequences profile 1:\n"); |
|---|
| 291 | for (i=0;i<nseqs;i++) |
|---|
| 292 | { |
|---|
| 293 | if (group[i+1] == 1) |
|---|
| 294 | { |
|---|
| 295 | if (debug>0) { |
|---|
| 296 | extern char **names; |
|---|
| 297 | fprintf(stdout,"%s\n",names[i+1]); |
|---|
| 298 | } |
|---|
| 299 | len = seqlen_array[i+1]; |
|---|
| 300 | alignment[nseqs1] = (char *) ckalloc( (len+2) * sizeof (char) ); |
|---|
| 301 | for (j=0;j<len;j++) |
|---|
| 302 | alignment[nseqs1][j] = seq_array[i+1][j+1]; |
|---|
| 303 | alignment[nseqs1][j] = ENDALN; |
|---|
| 304 | aln_len[nseqs1] = len; |
|---|
| 305 | aln_weight[nseqs1] = seq_weight[i]; |
|---|
| 306 | if (len > prf_length1) prf_length1 = len; |
|---|
| 307 | nseqs1++; |
|---|
| 308 | } |
|---|
| 309 | } |
|---|
| 310 | |
|---|
| 311 | /* |
|---|
| 312 | Make the second profile. |
|---|
| 313 | */ |
|---|
| 314 | prf_length2 = 0; |
|---|
| 315 | nseqs2 = 0; |
|---|
| 316 | if (debug>0) fprintf(stdout,"sequences profile 2:\n"); |
|---|
| 317 | for (i=0;i<nseqs;i++) |
|---|
| 318 | { |
|---|
| 319 | if (group[i+1] == 2) |
|---|
| 320 | { |
|---|
| 321 | if (debug>0) { |
|---|
| 322 | extern char **names; |
|---|
| 323 | fprintf(stdout,"%s\n",names[i+1]); |
|---|
| 324 | } |
|---|
| 325 | len = seqlen_array[i+1]; |
|---|
| 326 | alignment[nseqs1+nseqs2] = |
|---|
| 327 | (char *) ckalloc( (len+2) * sizeof (char) ); |
|---|
| 328 | for (j=0;j<len;j++) |
|---|
| 329 | alignment[nseqs1+nseqs2][j] = seq_array[i+1][j+1]; |
|---|
| 330 | alignment[nseqs1+nseqs2][j] = ENDALN; |
|---|
| 331 | aln_len[nseqs1+nseqs2] = len; |
|---|
| 332 | aln_weight[nseqs1+nseqs2] = seq_weight[i]; |
|---|
| 333 | if (len > prf_length2) prf_length2 = len; |
|---|
| 334 | nseqs2++; |
|---|
| 335 | } |
|---|
| 336 | } |
|---|
| 337 | |
|---|
| 338 | max_aln_length = prf_length1 + prf_length2+2; |
|---|
| 339 | |
|---|
| 340 | if (debug>0) fprintf(stdout,"average: %d\n",(pint)mat_avscore); |
|---|
| 341 | /* |
|---|
| 342 | calculate real length of profiles - removing end gaps! |
|---|
| 343 | */ |
|---|
| 344 | /* |
|---|
| 345 | is = ie = -1; |
|---|
| 346 | for (j=0; j<prf_length1; j++) |
|---|
| 347 | { |
|---|
| 348 | for (i=0;i<nseqs1;i++) |
|---|
| 349 | { |
|---|
| 350 | c = alignment[i][j]; |
|---|
| 351 | if ((c !=gap_pos1) && (c != gap_pos2)) |
|---|
| 352 | { |
|---|
| 353 | is = j; |
|---|
| 354 | break; |
|---|
| 355 | } |
|---|
| 356 | } |
|---|
| 357 | if (is != -1) break; |
|---|
| 358 | } |
|---|
| 359 | for (j=prf_length1-1; j>=0; j--) |
|---|
| 360 | { |
|---|
| 361 | for (i=0;i<nseqs1;i++) |
|---|
| 362 | { |
|---|
| 363 | c = alignment[i][j]; |
|---|
| 364 | if ((c !=gap_pos1) && (c != gap_pos2)) |
|---|
| 365 | { |
|---|
| 366 | ie = j; |
|---|
| 367 | break; |
|---|
| 368 | } |
|---|
| 369 | } |
|---|
| 370 | if (ie != -1) break; |
|---|
| 371 | } |
|---|
| 372 | len1 = ie-is; |
|---|
| 373 | |
|---|
| 374 | is = ie = -1; |
|---|
| 375 | for (j=0; j<prf_length2; j++) |
|---|
| 376 | { |
|---|
| 377 | for (i=nseqs1;i<nseqs1+nseqs2;i++) |
|---|
| 378 | { |
|---|
| 379 | c = alignment[i][j]; |
|---|
| 380 | if ((c !=gap_pos1) && (c != gap_pos2)) |
|---|
| 381 | { |
|---|
| 382 | is = j; |
|---|
| 383 | break; |
|---|
| 384 | } |
|---|
| 385 | } |
|---|
| 386 | if (is != -1) break; |
|---|
| 387 | } |
|---|
| 388 | for (j=prf_length2-1; j>=0; j--) |
|---|
| 389 | { |
|---|
| 390 | for (i=nseqs1;i<nseqs1+nseqs2;i++) |
|---|
| 391 | { |
|---|
| 392 | c = alignment[i][j]; |
|---|
| 393 | if ((c !=gap_pos1) && (c != gap_pos2)) |
|---|
| 394 | { |
|---|
| 395 | ie = j; |
|---|
| 396 | break; |
|---|
| 397 | } |
|---|
| 398 | } |
|---|
| 399 | if (ie != -1) break; |
|---|
| 400 | } |
|---|
| 401 | len2 = ie-is; |
|---|
| 402 | */ |
|---|
| 403 | len1=0; |
|---|
| 404 | for (i=0;i<nseqs1;i++) |
|---|
| 405 | { |
|---|
| 406 | is=0; |
|---|
| 407 | for (j=0; j<prf_length1; j++) |
|---|
| 408 | { |
|---|
| 409 | c = alignment[i][j]; |
|---|
| 410 | if ((c !=gap_pos1) && (c != gap_pos2)) is++; |
|---|
| 411 | } |
|---|
| 412 | if(is>len1) len1=is; |
|---|
| 413 | } |
|---|
| 414 | |
|---|
| 415 | len2=0; |
|---|
| 416 | for (i=0;i<nseqs2;i++) |
|---|
| 417 | { |
|---|
| 418 | is=0; |
|---|
| 419 | for (j=0; j<prf_length2; j++) |
|---|
| 420 | { |
|---|
| 421 | c = alignment[i][j]; |
|---|
| 422 | if ((c !=gap_pos1) && (c != gap_pos2)) is++; |
|---|
| 423 | } |
|---|
| 424 | if(is>len2) len2=is; |
|---|
| 425 | } |
|---|
| 426 | |
|---|
| 427 | if(len1==0 || len2==0) { |
|---|
| 428 | logmin=0; |
|---|
| 429 | logdiff=1.0; |
|---|
| 430 | } |
|---|
| 431 | else { |
|---|
| 432 | minlen = MIN(len1,len2); |
|---|
| 433 | if (logmin <= 1000) logmin = 0; |
|---|
| 434 | else |
|---|
| 435 | logmin = log((double)(minlen)); |
|---|
| 436 | if (len1<=len2) |
|---|
| 437 | logdiff = 1.0-log((double)((float)len1/(float)len2)); |
|---|
| 438 | else |
|---|
| 439 | logdiff = 1.0-log((double)((float)len2/(float)len1)); |
|---|
| 440 | } |
|---|
| 441 | /* |
|---|
| 442 | round logdiff to the nearest integer |
|---|
| 443 | */ |
|---|
| 444 | if ((logdiff-(int)logdiff) > 0.5) logdiff=ceil(logdiff); |
|---|
| 445 | else logdiff=floor(logdiff); |
|---|
| 446 | |
|---|
| 447 | if (debug>1) fprintf(stdout,"%d %d logmin %f logdiff %f\n", |
|---|
| 448 | (pint)len1,(pint)len2, logmin,logdiff); |
|---|
| 449 | |
|---|
| 450 | if (dnaflag) |
|---|
| 451 | { |
|---|
| 452 | gapcoef1 = gapcoef2 = 100.0 * gap_open *scale; |
|---|
| 453 | lencoef1 = lencoef2 = 100.0 * gap_extend *scale; |
|---|
| 454 | } |
|---|
| 455 | else |
|---|
| 456 | { |
|---|
| 457 | if (neg_matrix == TRUE) |
|---|
| 458 | { |
|---|
| 459 | if (mat_avscore <= 0) |
|---|
| 460 | gapcoef1 = gapcoef2 = 200.0 * (gap_open + logmin); |
|---|
| 461 | else |
|---|
| 462 | gapcoef1 = gapcoef2 = 200.0 * (gap_open + logmin); |
|---|
| 463 | lencoef1 = lencoef2 = 200.0 * gap_extend; |
|---|
| 464 | } |
|---|
| 465 | else |
|---|
| 466 | { |
|---|
| 467 | if (mat_avscore <= 0) |
|---|
| 468 | gapcoef1 = gapcoef2 = 100.0 * (float)(gap_open + logmin); |
|---|
| 469 | else |
|---|
| 470 | gapcoef1 = gapcoef2 = scale * mat_avscore * (float)(gap_open + logmin); |
|---|
| 471 | lencoef1 = lencoef2 = 100.0 * gap_extend; |
|---|
| 472 | } |
|---|
| 473 | if (len1>=len2) |
|---|
| 474 | { |
|---|
| 475 | gapcoef1 *= logdiff; |
|---|
| 476 | lencoef1 *= logdiff*15; |
|---|
| 477 | } |
|---|
| 478 | else |
|---|
| 479 | { |
|---|
| 480 | gapcoef2 *= logdiff; |
|---|
| 481 | lencoef2 *= logdiff*15; |
|---|
| 482 | } |
|---|
| 483 | } |
|---|
| 484 | |
|---|
| 485 | if (debug>0) |
|---|
| 486 | { |
|---|
| 487 | fprintf(stdout,"Gap Open1 %d Gap Open2 %d Gap Extend1 %d Gap Extend2 %d\n", |
|---|
| 488 | (pint)gapcoef1,(pint)gapcoef2, (pint)lencoef1,(pint)lencoef2); |
|---|
| 489 | fprintf(stdout,"Matrix %s\n", mtrxname); |
|---|
| 490 | } |
|---|
| 491 | |
|---|
| 492 | profile1 = (sint **) ckalloc( (prf_length1+2) * sizeof (sint *) ); |
|---|
| 493 | for(i=0; i<prf_length1+2; i++) |
|---|
| 494 | profile1[i] = (sint *) ckalloc( (LENCOL+2) * sizeof(sint) ); |
|---|
| 495 | |
|---|
| 496 | profile2 = (sint **) ckalloc( (prf_length2+2) * sizeof (sint *) ); |
|---|
| 497 | for(i=0; i<prf_length2+2; i++) |
|---|
| 498 | profile2[i] = (sint *) ckalloc( (LENCOL+2) * sizeof(sint) ); |
|---|
| 499 | |
|---|
| 500 | /* |
|---|
| 501 | calculate the Gap Coefficients. |
|---|
| 502 | */ |
|---|
| 503 | gaps = (sint *) ckalloc( (max_aln_length+1) * sizeof (sint) ); |
|---|
| 504 | |
|---|
| 505 | if (switch_profiles == FALSE) |
|---|
| 506 | calc_gap_coeff(alignment, gaps, profile1, (struct_penalties1 && use_ss1), gap_penalty_mask1, |
|---|
| 507 | (sint)0, nseqs1, prf_length1, gapcoef1, lencoef1); |
|---|
| 508 | else |
|---|
| 509 | calc_gap_coeff(alignment, gaps, profile1, (struct_penalties2 && use_ss2), gap_penalty_mask2, |
|---|
| 510 | (sint)0, nseqs1, prf_length1, gapcoef1, lencoef1); |
|---|
| 511 | /* |
|---|
| 512 | calculate the profile matrix. |
|---|
| 513 | */ |
|---|
| 514 | calc_prf1(profile1, alignment, gaps, matrix, |
|---|
| 515 | aln_weight, prf_length1, (sint)0, nseqs1); |
|---|
| 516 | |
|---|
| 517 | if (debug>4) |
|---|
| 518 | { |
|---|
| 519 | extern char *amino_acid_codes; |
|---|
| 520 | for (j=0;j<=max_aa;j++) |
|---|
| 521 | fprintf(stdout,"%c ", amino_acid_codes[j]); |
|---|
| 522 | fprintf(stdout,"\n"); |
|---|
| 523 | for (i=0;i<prf_length1;i++) |
|---|
| 524 | { |
|---|
| 525 | for (j=0;j<=max_aa;j++) |
|---|
| 526 | fprintf(stdout,"%d ", (pint)profile1[i+1][j]); |
|---|
| 527 | fprintf(stdout,"%d ", (pint)profile1[i+1][gap_pos1]); |
|---|
| 528 | fprintf(stdout,"%d ", (pint)profile1[i+1][gap_pos2]); |
|---|
| 529 | fprintf(stdout,"%d %d\n",(pint)profile1[i+1][GAPCOL],(pint)profile1[i+1][LENCOL]); |
|---|
| 530 | } |
|---|
| 531 | } |
|---|
| 532 | |
|---|
| 533 | /* |
|---|
| 534 | calculate the Gap Coefficients. |
|---|
| 535 | */ |
|---|
| 536 | |
|---|
| 537 | if (switch_profiles == FALSE) |
|---|
| 538 | calc_gap_coeff(alignment, gaps, profile2, (struct_penalties2 && use_ss2), gap_penalty_mask2, |
|---|
| 539 | nseqs1, nseqs1+nseqs2, prf_length2, gapcoef2, lencoef2); |
|---|
| 540 | else |
|---|
| 541 | calc_gap_coeff(alignment, gaps, profile2, (struct_penalties1 && use_ss1), gap_penalty_mask1, |
|---|
| 542 | nseqs1, nseqs1+nseqs2, prf_length2, gapcoef2, lencoef2); |
|---|
| 543 | /* |
|---|
| 544 | calculate the profile matrix. |
|---|
| 545 | */ |
|---|
| 546 | calc_prf2(profile2, alignment, aln_weight, |
|---|
| 547 | prf_length2, nseqs1, nseqs1+nseqs2); |
|---|
| 548 | |
|---|
| 549 | if (debug>4) |
|---|
| 550 | { |
|---|
| 551 | extern char *amino_acid_codes; |
|---|
| 552 | for (j=0;j<=max_aa;j++) |
|---|
| 553 | fprintf(stdout,"%c ", amino_acid_codes[j]); |
|---|
| 554 | fprintf(stdout,"\n"); |
|---|
| 555 | for (i=0;i<prf_length2;i++) |
|---|
| 556 | { |
|---|
| 557 | for (j=0;j<=max_aa;j++) |
|---|
| 558 | fprintf(stdout,"%d ", (pint)profile2[i+1][j]); |
|---|
| 559 | fprintf(stdout,"%d ", (pint)profile2[i+1][gap_pos1]); |
|---|
| 560 | fprintf(stdout,"%d ", (pint)profile2[i+1][gap_pos2]); |
|---|
| 561 | fprintf(stdout,"%d %d\n",(pint)profile2[i+1][GAPCOL],(pint)profile2[i+1][LENCOL]); |
|---|
| 562 | } |
|---|
| 563 | } |
|---|
| 564 | |
|---|
| 565 | aln_path1 = (char *) ckalloc( (max_aln_length+1) * sizeof(char) ); |
|---|
| 566 | aln_path2 = (char *) ckalloc( (max_aln_length+1) * sizeof(char) ); |
|---|
| 567 | |
|---|
| 568 | |
|---|
| 569 | /* |
|---|
| 570 | align the profiles |
|---|
| 571 | */ |
|---|
| 572 | /* use Myers and Miller to align two sequences */ |
|---|
| 573 | |
|---|
| 574 | last_print = 0; |
|---|
| 575 | print_ptr = 1; |
|---|
| 576 | |
|---|
| 577 | sb1 = sb2 = 0; |
|---|
| 578 | se1 = prf_length1; |
|---|
| 579 | se2 = prf_length2; |
|---|
| 580 | |
|---|
| 581 | HH = (lint *) ckalloc( (max_aln_length+1) * sizeof (lint) ); |
|---|
| 582 | DD = (lint *) ckalloc( (max_aln_length+1) * sizeof (lint) ); |
|---|
| 583 | RR = (lint *) ckalloc( (max_aln_length+1) * sizeof (lint) ); |
|---|
| 584 | SS = (lint *) ckalloc( (max_aln_length+1) * sizeof (lint) ); |
|---|
| 585 | gS = (lint *) ckalloc( (max_aln_length+1) * sizeof (lint) ); |
|---|
| 586 | displ = (sint *) ckalloc( (max_aln_length+1) * sizeof (sint) ); |
|---|
| 587 | |
|---|
| 588 | score = pdiff(sb1, sb2, se1-sb1, se2-sb2, 0, 0); |
|---|
| 589 | HH=ckfree((void *)HH); |
|---|
| 590 | DD=ckfree((void *)DD); |
|---|
| 591 | RR=ckfree((void *)RR); |
|---|
| 592 | SS=ckfree((void *)SS); |
|---|
| 593 | gS=ckfree((void *)gS); |
|---|
| 594 | |
|---|
| 595 | ptracepath( &alignment_len); |
|---|
| 596 | |
|---|
| 597 | displ=ckfree((void *)displ); |
|---|
| 598 | |
|---|
| 599 | add_ggaps(); |
|---|
| 600 | |
|---|
| 601 | for (i=0;i<prf_length1+2;i++) |
|---|
| 602 | profile1[i]=ckfree((void *)profile1[i]); |
|---|
| 603 | profile1=ckfree((void *)profile1); |
|---|
| 604 | |
|---|
| 605 | for (i=0;i<prf_length2+2;i++) |
|---|
| 606 | profile2[i]=ckfree((void *)profile2[i]); |
|---|
| 607 | profile2=ckfree((void *)profile2); |
|---|
| 608 | |
|---|
| 609 | prf_length1 = alignment_len; |
|---|
| 610 | |
|---|
| 611 | aln_path1=ckfree((void *)aln_path1); |
|---|
| 612 | aln_path2=ckfree((void *)aln_path2); |
|---|
| 613 | |
|---|
| 614 | NumSeq = 0; |
|---|
| 615 | for (j=0;j<nseqs;j++) |
|---|
| 616 | { |
|---|
| 617 | if (group[j+1] == 1) |
|---|
| 618 | { |
|---|
| 619 | seqlen_array[j+1] = prf_length1; |
|---|
| 620 | realloc_seq(j+1,prf_length1); |
|---|
| 621 | for (i=0;i<prf_length1;i++) |
|---|
| 622 | seq_array[j+1][i+1] = alignment[NumSeq][i]; |
|---|
| 623 | NumSeq++; |
|---|
| 624 | } |
|---|
| 625 | } |
|---|
| 626 | for (j=0;j<nseqs;j++) |
|---|
| 627 | { |
|---|
| 628 | if (group[j+1] == 2) |
|---|
| 629 | { |
|---|
| 630 | seqlen_array[j+1] = prf_length1; |
|---|
| 631 | seq_array[j+1] = (char *)realloc(seq_array[j+1], (prf_length1+2) * sizeof (char)); |
|---|
| 632 | realloc_seq(j+1,prf_length1); |
|---|
| 633 | for (i=0;i<prf_length1;i++) |
|---|
| 634 | seq_array[j+1][i+1] = alignment[NumSeq][i]; |
|---|
| 635 | NumSeq++; |
|---|
| 636 | } |
|---|
| 637 | } |
|---|
| 638 | |
|---|
| 639 | for (i=0;i<nseqs1+nseqs2;i++) |
|---|
| 640 | alignment[i]=ckfree((void *)alignment[i]); |
|---|
| 641 | alignment=ckfree((void *)alignment); |
|---|
| 642 | |
|---|
| 643 | aln_len=ckfree((void *)aln_len); |
|---|
| 644 | gaps=ckfree((void *)gaps); |
|---|
| 645 | |
|---|
| 646 | return(score/100); |
|---|
| 647 | } |
|---|
| 648 | |
|---|
| 649 | static void add_ggaps(void) |
|---|
| 650 | { |
|---|
| 651 | sint j; |
|---|
| 652 | sint i,ix; |
|---|
| 653 | sint len; |
|---|
| 654 | char *ta; |
|---|
| 655 | |
|---|
| 656 | ta = (char *) ckalloc( (alignment_len+1) * sizeof (char) ); |
|---|
| 657 | |
|---|
| 658 | for (j=0;j<nseqs1;j++) |
|---|
| 659 | { |
|---|
| 660 | ix = 0; |
|---|
| 661 | for (i=0;i<alignment_len;i++) |
|---|
| 662 | { |
|---|
| 663 | if (aln_path1[i] == 2) |
|---|
| 664 | { |
|---|
| 665 | if (ix < aln_len[j]) |
|---|
| 666 | ta[i] = alignment[j][ix]; |
|---|
| 667 | else |
|---|
| 668 | ta[i] = ENDALN; |
|---|
| 669 | ix++; |
|---|
| 670 | } |
|---|
| 671 | else if (aln_path1[i] == 1) |
|---|
| 672 | { |
|---|
| 673 | /* |
|---|
| 674 | insertion in first alignment... |
|---|
| 675 | */ |
|---|
| 676 | ta[i] = gap_pos1; |
|---|
| 677 | } |
|---|
| 678 | else |
|---|
| 679 | { |
|---|
| 680 | fprintf(stdout,"Error in aln_path\n"); |
|---|
| 681 | } |
|---|
| 682 | } |
|---|
| 683 | ta[i] = ENDALN; |
|---|
| 684 | |
|---|
| 685 | len = alignment_len; |
|---|
| 686 | alignment[j] = (char *)realloc(alignment[j], (len+2) * sizeof (char)); |
|---|
| 687 | for (i=0;i<len;i++) |
|---|
| 688 | alignment[j][i] = ta[i]; |
|---|
| 689 | alignment[j][i] = ENDALN; |
|---|
| 690 | aln_len[j] = len; |
|---|
| 691 | } |
|---|
| 692 | |
|---|
| 693 | for (j=nseqs1;j<nseqs1+nseqs2;j++) |
|---|
| 694 | { |
|---|
| 695 | ix = 0; |
|---|
| 696 | for (i=0;i<alignment_len;i++) |
|---|
| 697 | { |
|---|
| 698 | if (aln_path2[i] == 2) |
|---|
| 699 | { |
|---|
| 700 | if (ix < aln_len[j]) |
|---|
| 701 | ta[i] = alignment[j][ix]; |
|---|
| 702 | else |
|---|
| 703 | ta[i] = ENDALN; |
|---|
| 704 | ix++; |
|---|
| 705 | } |
|---|
| 706 | else if (aln_path2[i] == 1) |
|---|
| 707 | { |
|---|
| 708 | /* |
|---|
| 709 | insertion in second alignment... |
|---|
| 710 | */ |
|---|
| 711 | ta[i] = gap_pos1; |
|---|
| 712 | } |
|---|
| 713 | else |
|---|
| 714 | { |
|---|
| 715 | fprintf(stdout,"Error in aln_path\n"); |
|---|
| 716 | } |
|---|
| 717 | } |
|---|
| 718 | ta[i] = ENDALN; |
|---|
| 719 | |
|---|
| 720 | len = alignment_len; |
|---|
| 721 | alignment[j] = (char *) realloc(alignment[j], (len+2) * sizeof (char) ); |
|---|
| 722 | for (i=0;i<len;i++) |
|---|
| 723 | alignment[j][i] = ta[i]; |
|---|
| 724 | alignment[j][i] = ENDALN; |
|---|
| 725 | aln_len[j] = len; |
|---|
| 726 | } |
|---|
| 727 | |
|---|
| 728 | ta=ckfree((void *)ta); |
|---|
| 729 | |
|---|
| 730 | if (struct_penalties1 != NONE) |
|---|
| 731 | gap_penalty_mask1 = add_ggaps_mask(gap_penalty_mask1,alignment_len,aln_path1,aln_path2); |
|---|
| 732 | if (struct_penalties1 == SECST) |
|---|
| 733 | sec_struct_mask1 = add_ggaps_mask(sec_struct_mask1,alignment_len,aln_path1,aln_path2); |
|---|
| 734 | |
|---|
| 735 | if (struct_penalties2 != NONE) |
|---|
| 736 | gap_penalty_mask2 = add_ggaps_mask(gap_penalty_mask2,alignment_len,aln_path2,aln_path1); |
|---|
| 737 | if (struct_penalties2 == SECST) |
|---|
| 738 | sec_struct_mask2 = add_ggaps_mask(sec_struct_mask2,alignment_len,aln_path2,aln_path1); |
|---|
| 739 | |
|---|
| 740 | if (debug>0) |
|---|
| 741 | { |
|---|
| 742 | char c; |
|---|
| 743 | extern char *amino_acid_codes; |
|---|
| 744 | |
|---|
| 745 | for (i=0;i<nseqs1+nseqs2;i++) |
|---|
| 746 | { |
|---|
| 747 | for (j=0;j<alignment_len;j++) |
|---|
| 748 | { |
|---|
| 749 | if (alignment[i][j] == ENDALN) break; |
|---|
| 750 | else if ((alignment[i][j] == gap_pos1) || (alignment[i][j] == gap_pos2)) c = '-'; |
|---|
| 751 | else c = amino_acid_codes[alignment[i][j]]; |
|---|
| 752 | fprintf(stdout,"%c", c); |
|---|
| 753 | } |
|---|
| 754 | fprintf(stdout,"\n\n"); |
|---|
| 755 | } |
|---|
| 756 | } |
|---|
| 757 | |
|---|
| 758 | } |
|---|
| 759 | |
|---|
| 760 | static char * add_ggaps_mask(char *mask, int len, char *path1, char *path2) |
|---|
| 761 | { |
|---|
| 762 | int i,ix; |
|---|
| 763 | char *ta; |
|---|
| 764 | |
|---|
| 765 | ta = (char *) ckalloc( (len+1) * sizeof (char) ); |
|---|
| 766 | |
|---|
| 767 | ix = 0; |
|---|
| 768 | if (switch_profiles == FALSE) |
|---|
| 769 | { |
|---|
| 770 | for (i=0;i<len;i++) |
|---|
| 771 | { |
|---|
| 772 | if (path1[i] == 2) |
|---|
| 773 | { |
|---|
| 774 | ta[i] = mask[ix]; |
|---|
| 775 | ix++; |
|---|
| 776 | } |
|---|
| 777 | else if (path1[i] == 1) |
|---|
| 778 | ta[i] = gap_pos1; |
|---|
| 779 | } |
|---|
| 780 | } |
|---|
| 781 | else |
|---|
| 782 | { |
|---|
| 783 | for (i=0;i<len;i++) |
|---|
| 784 | { |
|---|
| 785 | if (path2[i] == 2) |
|---|
| 786 | { |
|---|
| 787 | ta[i] = mask[ix]; |
|---|
| 788 | ix++; |
|---|
| 789 | } |
|---|
| 790 | else if (path2[i] == 1) |
|---|
| 791 | ta[i] = gap_pos1; |
|---|
| 792 | } |
|---|
| 793 | } |
|---|
| 794 | mask = (char *)realloc(mask,(len+2) * sizeof (char)); |
|---|
| 795 | for (i=0;i<len;i++) |
|---|
| 796 | mask[i] = ta[i]; |
|---|
| 797 | mask[i] ='\0'; |
|---|
| 798 | |
|---|
| 799 | ta=ckfree((void *)ta); |
|---|
| 800 | |
|---|
| 801 | return(mask); |
|---|
| 802 | } |
|---|
| 803 | |
|---|
| 804 | static lint prfscore(sint n, sint m) |
|---|
| 805 | { |
|---|
| 806 | sint ix; |
|---|
| 807 | lint score; |
|---|
| 808 | |
|---|
| 809 | score = 0.0; |
|---|
| 810 | for (ix=0; ix<=max_aa; ix++) |
|---|
| 811 | { |
|---|
| 812 | score += (profile1[n][ix] * profile2[m][ix]); |
|---|
| 813 | } |
|---|
| 814 | score += (profile1[n][gap_pos1] * profile2[m][gap_pos1]); |
|---|
| 815 | score += (profile1[n][gap_pos2] * profile2[m][gap_pos2]); |
|---|
| 816 | return(score/10); |
|---|
| 817 | |
|---|
| 818 | } |
|---|
| 819 | |
|---|
| 820 | static void ptracepath(sint *alen) |
|---|
| 821 | { |
|---|
| 822 | sint i,j,k,pos,to_do; |
|---|
| 823 | |
|---|
| 824 | pos = 0; |
|---|
| 825 | |
|---|
| 826 | to_do=print_ptr-1; |
|---|
| 827 | |
|---|
| 828 | for(i=1;i<=to_do;++i) { |
|---|
| 829 | if (debug>1) fprintf(stdout,"%d ",(pint)displ[i]); |
|---|
| 830 | if(displ[i]==0) { |
|---|
| 831 | aln_path1[pos]=2; |
|---|
| 832 | aln_path2[pos]=2; |
|---|
| 833 | ++pos; |
|---|
| 834 | } |
|---|
| 835 | else { |
|---|
| 836 | if((k=displ[i])>0) { |
|---|
| 837 | for(j=0;j<=k-1;++j) { |
|---|
| 838 | aln_path2[pos+j]=2; |
|---|
| 839 | aln_path1[pos+j]=1; |
|---|
| 840 | } |
|---|
| 841 | pos += k; |
|---|
| 842 | } |
|---|
| 843 | else { |
|---|
| 844 | k = (displ[i]<0) ? displ[i] * -1 : displ[i]; |
|---|
| 845 | for(j=0;j<=k-1;++j) { |
|---|
| 846 | aln_path1[pos+j]=2; |
|---|
| 847 | aln_path2[pos+j]=1; |
|---|
| 848 | } |
|---|
| 849 | pos += k; |
|---|
| 850 | } |
|---|
| 851 | } |
|---|
| 852 | } |
|---|
| 853 | if (debug>1) fprintf(stdout,"\n"); |
|---|
| 854 | |
|---|
| 855 | (*alen) = pos; |
|---|
| 856 | |
|---|
| 857 | } |
|---|
| 858 | |
|---|
| 859 | static void pdel(sint k) |
|---|
| 860 | { |
|---|
| 861 | if(last_print<0) |
|---|
| 862 | last_print = displ[print_ptr-1] -= k; |
|---|
| 863 | else |
|---|
| 864 | last_print = displ[print_ptr++] = -(k); |
|---|
| 865 | } |
|---|
| 866 | |
|---|
| 867 | static void padd(sint k) |
|---|
| 868 | { |
|---|
| 869 | |
|---|
| 870 | if(last_print<0) { |
|---|
| 871 | displ[print_ptr-1] = k; |
|---|
| 872 | displ[print_ptr++] = last_print; |
|---|
| 873 | } |
|---|
| 874 | else |
|---|
| 875 | last_print = displ[print_ptr++] = k; |
|---|
| 876 | } |
|---|
| 877 | |
|---|
| 878 | static void palign(void) |
|---|
| 879 | { |
|---|
| 880 | displ[print_ptr++] = last_print = 0; |
|---|
| 881 | } |
|---|
| 882 | |
|---|
| 883 | |
|---|
| 884 | static lint pdiff(sint A,sint B,sint M,sint N,sint go1, sint go2) |
|---|
| 885 | { |
|---|
| 886 | sint midi,midj,type; |
|---|
| 887 | lint midh; |
|---|
| 888 | |
|---|
| 889 | static lint t, tl, g, h; |
|---|
| 890 | |
|---|
| 891 | { static sint i,j; |
|---|
| 892 | static lint hh, f, e, s; |
|---|
| 893 | |
|---|
| 894 | /* Boundary cases: M <= 1 or N == 0 */ |
|---|
| 895 | if (debug>2) fprintf(stdout,"A %d B %d M %d N %d midi %d go1 %d go2 %d\n", |
|---|
| 896 | (pint)A,(pint)B,(pint)M,(pint)N,(pint)M/2,(pint)go1,(pint)go2); |
|---|
| 897 | |
|---|
| 898 | /* if sequence B is empty.... */ |
|---|
| 899 | |
|---|
| 900 | if(N<=0) { |
|---|
| 901 | |
|---|
| 902 | /* if sequence A is not empty.... */ |
|---|
| 903 | |
|---|
| 904 | if(M>0) { |
|---|
| 905 | |
|---|
| 906 | /* delete residues A[1] to A[M] */ |
|---|
| 907 | |
|---|
| 908 | pdel(M); |
|---|
| 909 | } |
|---|
| 910 | return(-gap_penalty1(A,B,M)); |
|---|
| 911 | } |
|---|
| 912 | |
|---|
| 913 | /* if sequence A is empty.... */ |
|---|
| 914 | |
|---|
| 915 | if(M<=1) { |
|---|
| 916 | if(M<=0) { |
|---|
| 917 | |
|---|
| 918 | /* insert residues B[1] to B[N] */ |
|---|
| 919 | |
|---|
| 920 | padd(N); |
|---|
| 921 | return(-gap_penalty2(A,B,N)); |
|---|
| 922 | } |
|---|
| 923 | |
|---|
| 924 | /* if sequence A has just one residue.... */ |
|---|
| 925 | |
|---|
| 926 | if (go1 == 0) |
|---|
| 927 | midh = -gap_penalty1(A+1,B+1,N); |
|---|
| 928 | else |
|---|
| 929 | midh = -gap_penalty2(A+1,B,1)-gap_penalty1(A+1,B+1,N); |
|---|
| 930 | if (go2 == 0) |
|---|
| 931 | hh = -gap_penalty1(A,B+1,N); |
|---|
| 932 | else |
|---|
| 933 | hh = -gap_penalty1(A,B+1,N)-gap_penalty2(A+1,B+N,1); |
|---|
| 934 | if(hh>midh) midh = hh; |
|---|
| 935 | midj = 0; |
|---|
| 936 | for(j=1;j<=N;j++) { |
|---|
| 937 | hh = -gap_penalty1(A,B+1,j-1) + prfscore(A+1,B+j) |
|---|
| 938 | -gap_penalty1(A+1,B+j+1,N-j); |
|---|
| 939 | if(hh>midh) { |
|---|
| 940 | midh = hh; |
|---|
| 941 | midj = j; |
|---|
| 942 | } |
|---|
| 943 | } |
|---|
| 944 | |
|---|
| 945 | if(midj==0) { |
|---|
| 946 | pdel(1); |
|---|
| 947 | padd(N); |
|---|
| 948 | } |
|---|
| 949 | else { |
|---|
| 950 | if(midj>1) padd(midj-1); |
|---|
| 951 | palign(); |
|---|
| 952 | if(midj<N) padd(N-midj); |
|---|
| 953 | } |
|---|
| 954 | return midh; |
|---|
| 955 | } |
|---|
| 956 | |
|---|
| 957 | |
|---|
| 958 | /* Divide sequence A in half: midi */ |
|---|
| 959 | |
|---|
| 960 | midi = M / 2; |
|---|
| 961 | |
|---|
| 962 | /* In a forward phase, calculate all HH[j] and HH[j] */ |
|---|
| 963 | |
|---|
| 964 | HH[0] = 0.0; |
|---|
| 965 | t = -open_penalty1(A,B+1); |
|---|
| 966 | tl = -ext_penalty1(A,B+1); |
|---|
| 967 | for(j=1;j<=N;j++) { |
|---|
| 968 | HH[j] = t = t+tl; |
|---|
| 969 | DD[j] = t-open_penalty2(A+1,B+j); |
|---|
| 970 | } |
|---|
| 971 | |
|---|
| 972 | if (go1 == 0) t = 0; |
|---|
| 973 | else t = -open_penalty2(A+1,B); |
|---|
| 974 | tl = -ext_penalty2(A+1,B); |
|---|
| 975 | for(i=1;i<=midi;i++) { |
|---|
| 976 | s = HH[0]; |
|---|
| 977 | HH[0] = hh = t = t+tl; |
|---|
| 978 | f = t-open_penalty1(A+i,B+1); |
|---|
| 979 | |
|---|
| 980 | for(j=1;j<=N;j++) { |
|---|
| 981 | g = open_penalty1(A+i,B+j); |
|---|
| 982 | h = ext_penalty1(A+i,B+j); |
|---|
| 983 | if ((hh=hh-g-h) > (f=f-h)) f=hh; |
|---|
| 984 | g = open_penalty2(A+i,B+j); |
|---|
| 985 | h = ext_penalty2(A+i,B+j); |
|---|
| 986 | if ((hh=HH[j]-g-h) > (e=DD[j]-h)) e=hh; |
|---|
| 987 | hh = s + prfscore(A+i, B+j); |
|---|
| 988 | if (f>hh) hh = f; |
|---|
| 989 | if (e>hh) hh = e; |
|---|
| 990 | |
|---|
| 991 | s = HH[j]; |
|---|
| 992 | HH[j] = hh; |
|---|
| 993 | DD[j] = e; |
|---|
| 994 | |
|---|
| 995 | } |
|---|
| 996 | } |
|---|
| 997 | |
|---|
| 998 | DD[0]=HH[0]; |
|---|
| 999 | |
|---|
| 1000 | /* In a reverse phase, calculate all RR[j] and SS[j] */ |
|---|
| 1001 | |
|---|
| 1002 | RR[N]=0.0; |
|---|
| 1003 | tl = 0.0; |
|---|
| 1004 | for(j=N-1;j>=0;j--) { |
|---|
| 1005 | g = -open_penalty1(A+M,B+j+1); |
|---|
| 1006 | tl -= ext_penalty1(A+M,B+j+1); |
|---|
| 1007 | RR[j] = g+tl; |
|---|
| 1008 | SS[j] = RR[j]-open_penalty2(A+M,B+j); |
|---|
| 1009 | gS[j] = open_penalty2(A+M,B+j); |
|---|
| 1010 | } |
|---|
| 1011 | |
|---|
| 1012 | tl = 0.0; |
|---|
| 1013 | for(i=M-1;i>=midi;i--) { |
|---|
| 1014 | s = RR[N]; |
|---|
| 1015 | if (go2 == 0) g = 0; |
|---|
| 1016 | else g = -open_penalty2(A+i+1,B+N); |
|---|
| 1017 | tl -= ext_penalty2(A+i+1,B+N); |
|---|
| 1018 | RR[N] = hh = g+tl; |
|---|
| 1019 | t = open_penalty1(A+i,B+N); |
|---|
| 1020 | f = RR[N]-t; |
|---|
| 1021 | |
|---|
| 1022 | for(j=N-1;j>=0;j--) { |
|---|
| 1023 | g = open_penalty1(A+i,B+j+1); |
|---|
| 1024 | h = ext_penalty1(A+i,B+j+1); |
|---|
| 1025 | if ((hh=hh-g-h) > (f=f-h-g+t)) f=hh; |
|---|
| 1026 | t = g; |
|---|
| 1027 | g = open_penalty2(A+i+1,B+j); |
|---|
| 1028 | h = ext_penalty2(A+i+1,B+j); |
|---|
| 1029 | hh=RR[j]-g-h; |
|---|
| 1030 | if (i==(M-1)) { |
|---|
| 1031 | e=SS[j]-h; |
|---|
| 1032 | } |
|---|
| 1033 | else { |
|---|
| 1034 | e=SS[j]-h-g+open_penalty2(A+i+2,B+j); |
|---|
| 1035 | gS[j] = g; |
|---|
| 1036 | } |
|---|
| 1037 | if (hh > e) e=hh; |
|---|
| 1038 | hh = s + prfscore(A+i+1, B+j+1); |
|---|
| 1039 | if (f>hh) hh = f; |
|---|
| 1040 | if (e>hh) hh = e; |
|---|
| 1041 | |
|---|
| 1042 | s = RR[j]; |
|---|
| 1043 | RR[j] = hh; |
|---|
| 1044 | SS[j] = e; |
|---|
| 1045 | |
|---|
| 1046 | } |
|---|
| 1047 | } |
|---|
| 1048 | SS[N]=RR[N]; |
|---|
| 1049 | if (go2 == 0) gS[N] = 0; |
|---|
| 1050 | else gS[N] = open_penalty2(A+midi+1,B+N); |
|---|
| 1051 | |
|---|
| 1052 | /* find midj, such that HH[j]+RR[j] or DD[j]+SS[j]+gap is the maximum */ |
|---|
| 1053 | |
|---|
| 1054 | midh=HH[0]+RR[0]; |
|---|
| 1055 | midj=0; |
|---|
| 1056 | type=1; |
|---|
| 1057 | for(j=0;j<=N;j++) { |
|---|
| 1058 | hh = HH[j] + RR[j]; |
|---|
| 1059 | if(hh>=midh) |
|---|
| 1060 | if(hh>midh || (HH[j]!=DD[j] && RR[j]==SS[j])) { |
|---|
| 1061 | midh=hh; |
|---|
| 1062 | midj=j; |
|---|
| 1063 | } |
|---|
| 1064 | } |
|---|
| 1065 | |
|---|
| 1066 | for(j=N;j>=0;j--) { |
|---|
| 1067 | g = open_penalty2(A+midi+1,B+j); |
|---|
| 1068 | hh = DD[j] + SS[j] + gS[j]; |
|---|
| 1069 | if(hh>midh) { |
|---|
| 1070 | midh=hh; |
|---|
| 1071 | midj=j; |
|---|
| 1072 | type=2; |
|---|
| 1073 | } |
|---|
| 1074 | } |
|---|
| 1075 | } |
|---|
| 1076 | |
|---|
| 1077 | /* Conquer recursively around midpoint */ |
|---|
| 1078 | |
|---|
| 1079 | |
|---|
| 1080 | if(type==1) { /* Type 1 gaps */ |
|---|
| 1081 | if (debug>2) fprintf(stdout,"Type 1,1: midj %d\n",(pint)midj); |
|---|
| 1082 | pdiff(A,B,midi,midj,1,1); |
|---|
| 1083 | if (debug>2) fprintf(stdout,"Type 1,2: midj %d\n",(pint)midj); |
|---|
| 1084 | pdiff(A+midi,B+midj,M-midi,N-midj,1,1); |
|---|
| 1085 | } |
|---|
| 1086 | else { |
|---|
| 1087 | if (debug>2) fprintf(stdout,"Type 2,1: midj %d\n",(pint)midj); |
|---|
| 1088 | pdiff(A,B,midi-1,midj,1, 0); |
|---|
| 1089 | pdel(2); |
|---|
| 1090 | if (debug>2) fprintf(stdout,"Type 2,2: midj %d\n",(pint)midj); |
|---|
| 1091 | pdiff(A+midi+1,B+midj,M-midi-1,N-midj,0,1); |
|---|
| 1092 | } |
|---|
| 1093 | |
|---|
| 1094 | return midh; /* Return the score of the best alignment */ |
|---|
| 1095 | } |
|---|
| 1096 | |
|---|
| 1097 | /* calculate the score for opening a gap at residues A[i] and B[j] */ |
|---|
| 1098 | |
|---|
| 1099 | static sint open_penalty1(sint i, sint j) |
|---|
| 1100 | { |
|---|
| 1101 | sint g; |
|---|
| 1102 | |
|---|
| 1103 | if (i==0 || i==prf_length1) return(0); |
|---|
| 1104 | |
|---|
| 1105 | g = profile2[j][GAPCOL] + profile1[i][GAPCOL]; |
|---|
| 1106 | return(g); |
|---|
| 1107 | } |
|---|
| 1108 | |
|---|
| 1109 | /* calculate the score for extending an existing gap at A[i] and B[j] */ |
|---|
| 1110 | |
|---|
| 1111 | static sint ext_penalty1(sint i, sint j) |
|---|
| 1112 | { |
|---|
| 1113 | sint h; |
|---|
| 1114 | |
|---|
| 1115 | if (i==0 || i==prf_length1) return(0); |
|---|
| 1116 | |
|---|
| 1117 | h = profile2[j][LENCOL]; |
|---|
| 1118 | return(h); |
|---|
| 1119 | } |
|---|
| 1120 | |
|---|
| 1121 | /* calculate the score for a gap of length k, at residues A[i] and B[j] */ |
|---|
| 1122 | |
|---|
| 1123 | static sint gap_penalty1(sint i, sint j, sint k) |
|---|
| 1124 | { |
|---|
| 1125 | sint ix; |
|---|
| 1126 | sint gp; |
|---|
| 1127 | sint g, h = 0; |
|---|
| 1128 | |
|---|
| 1129 | if (k <= 0) return(0); |
|---|
| 1130 | if (i==0 || i==prf_length1) return(0); |
|---|
| 1131 | |
|---|
| 1132 | g = profile2[j][GAPCOL] + profile1[i][GAPCOL]; |
|---|
| 1133 | for (ix=0;ix<k && ix+j<prf_length2;ix++) |
|---|
| 1134 | h = profile2[ix+j][LENCOL]; |
|---|
| 1135 | |
|---|
| 1136 | gp = g + h * k; |
|---|
| 1137 | return(gp); |
|---|
| 1138 | } |
|---|
| 1139 | /* calculate the score for opening a gap at residues A[i] and B[j] */ |
|---|
| 1140 | |
|---|
| 1141 | static sint open_penalty2(sint i, sint j) |
|---|
| 1142 | { |
|---|
| 1143 | sint g; |
|---|
| 1144 | |
|---|
| 1145 | if (j==0 || j==prf_length2) return(0); |
|---|
| 1146 | |
|---|
| 1147 | g = profile1[i][GAPCOL] + profile2[j][GAPCOL]; |
|---|
| 1148 | return(g); |
|---|
| 1149 | } |
|---|
| 1150 | |
|---|
| 1151 | /* calculate the score for extending an existing gap at A[i] and B[j] */ |
|---|
| 1152 | |
|---|
| 1153 | static sint ext_penalty2(sint i, sint j) |
|---|
| 1154 | { |
|---|
| 1155 | sint h; |
|---|
| 1156 | |
|---|
| 1157 | if (j==0 || j==prf_length2) return(0); |
|---|
| 1158 | |
|---|
| 1159 | h = profile1[i][LENCOL]; |
|---|
| 1160 | return(h); |
|---|
| 1161 | } |
|---|
| 1162 | |
|---|
| 1163 | /* calculate the score for a gap of length k, at residues A[i] and B[j] */ |
|---|
| 1164 | |
|---|
| 1165 | static sint gap_penalty2(sint i, sint j, sint k) |
|---|
| 1166 | { |
|---|
| 1167 | sint ix; |
|---|
| 1168 | sint gp; |
|---|
| 1169 | sint g, h = 0; |
|---|
| 1170 | |
|---|
| 1171 | if (k <= 0) return(0); |
|---|
| 1172 | if (j==0 || j==prf_length2) return(0); |
|---|
| 1173 | |
|---|
| 1174 | g = profile1[i][GAPCOL] + profile2[j][GAPCOL]; |
|---|
| 1175 | for (ix=0;ix<k && ix+i<prf_length1;ix++) |
|---|
| 1176 | h = profile1[ix+i][LENCOL]; |
|---|
| 1177 | |
|---|
| 1178 | gp = g + h * k; |
|---|
| 1179 | return(gp); |
|---|
| 1180 | } |
|---|
