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