annotate clustalomega/clustal-omega-1.0.2/src/clustal/ktuple_pair.c @ 1:bc707542e5de

Uploaded
author clustalomega
date Thu, 21 Jul 2011 13:35:08 -0400
parents
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
1
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
1 /* -*- mode: c; tab-width: 4; c-basic-offset: 4; indent-tabs-mode: nil -*- */
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
2
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
3 /*********************************************************************
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
4 * Clustal Omega - Multiple sequence alignment
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
5 *
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
6 * Copyright (C) 2010 University College Dublin
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
7 *
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
8 * Clustal-Omega is free software; you can redistribute it and/or
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
9 * modify it under the terms of the GNU General Public License as
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
10 * published by the Free Software Foundation; either version 2 of the
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
11 * License, or (at your option) any later version.
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
12 *
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
13 * This file is part of Clustal-Omega.
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
14 *
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
15 ********************************************************************/
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
16
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
17 /*
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
18 * RCS $Id: ktuple_pair.c 230 2011-04-09 15:37:50Z andreas $
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
19 *
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
20 *
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
21 * K-Tuple code for pairwise alignment (Wilbur and Lipman, 1983; PMID
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
22 * 6572363). Most code taken from showpair.c (Clustal 1.83)
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
23 * DD: some functions now have lots of parameters as static variables
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
24 * were removed to make code OpenMP-friendly
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
25 *
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
26 */
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
27
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
28 #ifdef HAVE_CONFIG_H
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
29 #include "config.h"
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
30 #endif
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
31
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
32
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
33 #include <stdio.h>
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
34 #include <string.h>
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
35 #include <ctype.h>
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
36 #include <stdlib.h>
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
37 #include <math.h>
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
38 #include <assert.h>
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
39
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
40 #ifdef HAVE_OPENMP
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
41 #include <omp.h>
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
42 #endif
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
43
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
44 #include "squid/squid.h"
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
45 #include "util.h"
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
46 #include "symmatrix.h"
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
47 #include "ktuple_pair.h"
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
48 #include "log.h"
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
49 #include "progress.h"
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
50
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
51 #define END_MARK -3 /* see interface.c in 1.83 */
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
52 #define NUMRES 32 /* max size of comparison matrix */
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
53
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
54 /* see notes below */
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
55 #undef SORT_LAST_ELEMENT_AS_WELL
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
56
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
57 /* gap_pos1 = NUMRES-2; /@ code for gaps inserted by clustalw @/ */
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
58 static const int GAP_POS2 = NUMRES-1; /* code for gaps already in alignment */
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
59 static bool DNAFLAG = FALSE;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
60
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
61 static const char *AMINO_ACID_CODES = "ABCDEFGHIKLMNPQRSTUVWXYZ-";
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
62 static const char *NUCLEIC_ACID_CODES = "ACGTUN-";
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
63 /* As far as I understand the gap symbol should not be necessary here,
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
64 * because we use isgap for testing later anyway. But changing this,
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
65 * will affect max_res_code and max_nuc as well. So I leave it for now
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
66 * as it is. AW
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
67 */
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
68
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
69 static bool percent = TRUE;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
70
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
71 static void make_ptrs(int *tptr, int *pl, const int naseq, const int l, const int ktup, const int max_res_code, char **seq_array);
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
72 static void put_frag(const int fs, const int v1, const int v2, const int flen, const int curr_frag, int *next, int *maxsf, int **accum);
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
73 static bool frag_rel_pos(int a1, int b1, int a2, int b2, int ktup);
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
74 static void des_quick_sort(int *array1, int *array2, const int array_size);
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
75 static void pair_align(int seq_no, int l1, int l2, int max_res_code, ktuple_param_t *aln_param,
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
76 char **seq_array, int *maxsf, int **accum, int max_aln_length,
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
77 int *zza, int *zzb, int *zzc, int *zzd);
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
78 static void encode(char *seq, char *naseq, int l, const char *res_codes);
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
79 static int res_index(const char *lookup, char c);
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
80
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
81
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
82 typedef struct {
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
83 int i1;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
84 int i2;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
85 } two_ints_t;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
86
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
87
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
88
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
89 /* default ktuple pairwise alignment parameters
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
90 *
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
91 */
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
92 /* protein
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
93 */
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
94 /* designated initializer */
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
95 const ktuple_param_t default_protein_param = {
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
96 .ktup = 1,
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
97 .wind_gap = 3,
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
98 .signif = 5,
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
99 .window = 5,
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
100 };
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
101 /* dna
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
102 */
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
103 /* designated initializer */
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
104 const ktuple_param_t default_dna_param = {
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
105 .ktup = 2,
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
106 .wind_gap = 5,
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
107 .signif = 4,
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
108 .window = 4,
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
109 };
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
110
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
111
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
112 /**
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
113 * note: naseq should be unit-offset
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
114 */
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
115 static void
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
116 encode(char *seq, char *naseq, int l, const char *res_codes)
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
117 {
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
118 /* code seq as ints .. use GAP_POS2 for gap */
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
119 register int i;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
120 bool seq_contains_unknown_char = FALSE;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
121 /*LOG_DEBUG("seq=%s naseq=%p l=%d", &(seq[1]), naseq, l); */
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
122
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
123
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
124 for (i=1; i<=l; i++) {
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
125 char res = toupper(seq[i]);
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
126 if (isgap(res)) {
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
127 naseq[i] = GAP_POS2; /* gap in input */
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
128 } else {
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
129 naseq[i] = res_index(res_codes, res);
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
130 }
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
131
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
132 /*LOG_DEBUG("Character '%c' at pos %d", res, i);*/
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
133 if (-1 == naseq[i]) {
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
134 seq_contains_unknown_char = TRUE;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
135 /*LOG_DEBUG("Unknown character '%c' at pos %d", res, i);*/
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
136 }
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
137 /*LOG_DEBUG("na_seq[%d]=%d", i, naseq[i]);*/
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
138 }
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
139
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
140 if (TRUE == seq_contains_unknown_char)
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
141 Log(&rLog, LOG_WARN, "Unknown character in seq '%s'", &(seq[1]));
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
142
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
143 naseq[i] = END_MARK;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
144
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
145 return;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
146 }
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
147 /* end of encode */
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
148
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
149
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
150 /**
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
151 *
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
152 */
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
153 static int
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
154 res_index(const char *t, char c)
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
155 {
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
156 register int i;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
157 for (i=0; t[i] && t[i] != c; i++)
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
158 ;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
159 if (t[i]) {
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
160 return (i);
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
161 } else {
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
162 return -1;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
163 }
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
164 }
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
165 /* end of res_index */
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
166
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
167
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
168 /**
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
169 *
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
170 */
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
171 static void
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
172 make_ptrs(int *tptr, int *pl, const int naseq, const int l, const int ktup, const int max_res_code, char **seq_array)
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
173 {
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
174 /* FIXME make 10 a constant and give it a nice name */
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
175 static int a[10];
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
176 int i, j, code, flag;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
177 char residue;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
178 const int limit = (int) pow((double)(max_res_code+1),(double)ktup);
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
179
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
180 for (i=1;i<=ktup;i++)
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
181 a[i] = (int) pow((double)(max_res_code+1),(double)(i-1));
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
182
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
183 for (i=1; i<=limit; ++i)
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
184 pl[i]=0;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
185 for (i=1; i<=l; ++i)
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
186 tptr[i]=0;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
187
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
188 for (i=1; i<=(l-ktup+1); ++i) {
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
189 code=0;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
190 flag=FALSE;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
191 for (j=1; j<=ktup; ++j) {
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
192 /* Log(&rLog, LOG_FORCED_DEBUG, "naseq=%d i=%d j=%d seq_array[naseq]=%p",
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
193 * naseq, i, j, seq_array[naseq]);
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
194 */
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
195 residue = seq_array[naseq][i+j-1];
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
196 /* Log(&rLog, LOG_FORCED_DEBUG, "residue = %d", residue); */
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
197 if ((residue<0) || (residue > max_res_code)){
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
198 flag=TRUE;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
199 break;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
200 }
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
201 code += ((residue) * a[j]);
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
202 }
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
203 if (flag)
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
204 continue;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
205 ++code;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
206 if (0 != pl[code])
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
207 tptr[i] =pl[code];
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
208 pl[code] = i;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
209 }
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
210
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
211 return;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
212 }
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
213 /* end of make_ptrs */
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
214
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
215
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
216 /**
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
217 *
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
218 * FIXME Why hardcoding of 5?
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
219 */
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
220 static void
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
221 put_frag(const int fs, const int v1, const int v2, const int flen, const int curr_frag, int *next, int *maxsf, int **accum)
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
222 {
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
223 int end;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
224 accum[0][curr_frag]=fs;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
225 accum[1][curr_frag]=v1;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
226 accum[2][curr_frag]=v2;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
227 accum[3][curr_frag]=flen;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
228
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
229 if (!*maxsf) {
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
230 *maxsf=1;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
231 accum[4][curr_frag]=0;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
232 return;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
233 }
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
234
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
235 if (fs >= accum[0][*maxsf]) {
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
236 accum[4][curr_frag]=*maxsf;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
237 *maxsf=curr_frag;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
238 return;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
239 } else {
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
240 *next=*maxsf;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
241 while (TRUE) {
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
242 end=*next;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
243 *next=accum[4][*next];
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
244 if (fs>=accum[0][*next])
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
245 break;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
246 }
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
247 accum[4][curr_frag]=*next;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
248 accum[4][end]=curr_frag;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
249 }
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
250
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
251 return;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
252 }
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
253 /* end of put_frag */
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
254
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
255
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
256 /**
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
257 *
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
258 */
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
259 static bool
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
260 frag_rel_pos(int a1, int b1, int a2, int b2, int ktup)
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
261 {
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
262 if (a1-b1 == a2-b2) {
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
263 if (a2<a1) {
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
264 return TRUE;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
265 }
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
266 } else {
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
267 if (a2+ktup-1<a1 && b2+ktup-1<b1) {
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
268 return TRUE;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
269 }
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
270 }
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
271 return FALSE;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
272 }
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
273 /* end of frag_rel_pos */
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
274
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
275
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
276
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
277
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
278 /**
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
279 *
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
280 * @note: This is together with des_quick_sort most time consuming
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
281 * routine according to gprof on r110. Tried to replace it with qsort
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
282 * and/or QSortAndTrackIndex(), which is always slower! So we keep the
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
283 * original.
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
284 *
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
285 * Original doc: Quicksort routine, adapted from chapter 4, page 115
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
286 * of software tools by Kernighan and Plauger, (1986). Sort the
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
287 * elements of array1 and sort the elements of array2 accordingly
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
288 *
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
289 * There might be a bug here. The original function apparently never
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
290 * touches the last element and keeps it as is. Tried to fix this (see
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
291 * SORT_LAST_ELEMENT_AS_WELL) which gives slightly worse performance
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
292 * (-0.5% on BB). My fix might not be working or it's not a bug at
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
293 * all...
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
294 *
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
295 *
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
296 *
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
297 */
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
298 static void
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
299 des_quick_sort(int *array1, int *array2, const int array_size)
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
300 {
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
301 int temp1, temp2;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
302 int p, pivlin;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
303 int i, j;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
304 int lst[50], ust[50]; /* the maximum no. of elements must be*/
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
305 /* < log(base2) of 50 */
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
306
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
307 #if 0
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
308 for (i=1; i<=array_size; i++) {
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
309 Log(&rLog, LOG_FORCED_DEBUG, "b4 sort array1[%d]=%d array2[%d]=%d", i, array1[i], i, array2[i]);
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
310 }
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
311 #endif
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
312 lst[1] = 1;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
313
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
314 #ifdef SORT_LAST_ELEMENT_AS_WELL
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
315 ust[1] = array_size;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
316 #else
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
317 /* original */
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
318 ust[1] = array_size-1;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
319 #endif
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
320 p = 1;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
321
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
322
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
323 while (p > 0) {
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
324 if (lst[p] >= ust[p]) {
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
325 p--;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
326 } else {
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
327 i = lst[p] - 1;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
328 j = ust[p];
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
329 pivlin = array1[j];
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
330 while (i < j) {
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
331 for (i=i+1; array1[i] < pivlin; i++)
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
332 ;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
333 for (j=j-1; j > i; j--)
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
334 if (array1[j] <= pivlin) break;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
335 if (i < j) {
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
336 temp1 = array1[i];
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
337 array1[i] = array1[j];
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
338 array1[j] = temp1;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
339
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
340 temp2 = array2[i];
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
341 array2[i] = array2[j];
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
342 array2[j] = temp2;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
343 }
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
344 }
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
345
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
346 j = ust[p];
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
347
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
348 temp1 = array1[i];
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
349 array1[i] = array1[j];
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
350 array1[j] = temp1;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
351
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
352 temp2 = array2[i];
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
353 array2[i] = array2[j];
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
354 array2[j] = temp2;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
355
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
356 if (i-lst[p] < ust[p] - i) {
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
357 lst[p+1] = lst[p];
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
358 ust[p+1] = i - 1;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
359 lst[p] = i + 1;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
360
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
361 } else {
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
362 lst[p+1] = i + 1;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
363 ust[p+1] = ust[p];
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
364 ust[p] = i - 1;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
365 }
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
366 p = p + 1;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
367 }
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
368 }
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
369
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
370 #if 0
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
371 for (i=1; i<=array_size; i++) {
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
372 Log(&rLog, LOG_FORCED_DEBUG, "after sort array1[%d]=%d array2[%d]=%d", i, array1[i], i, array2[i]);
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
373 }
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
374 #endif
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
375
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
376 return;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
377 }
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
378 /* end of des_quick_sort */
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
379
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
380
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
381
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
382 /**
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
383 *
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
384 * FIXME together with des_quick_sort most time consuming routine
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
385 * according to gprof on r110
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
386 *
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
387 */
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
388 static void
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
389 pair_align(int seq_no, int l1, int l2, int max_res_code, ktuple_param_t *aln_param,
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
390 char **seq_array, int *maxsf, int **accum, int max_aln_length,
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
391 int *zza, int *zzb, int *zzc, int *zzd)
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
392 {
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
393 int next; /* forrmerly static */
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
394 int pot[8],i, j, l, m, flag, limit, pos, vn1, vn2, flen, osptr, fs;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
395 int tv1, tv2, encrypt, subt1, subt2, rmndr;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
396 char residue;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
397 int *diag_index;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
398 int *displ;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
399 char *slopes;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
400 int curr_frag;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
401 const int tl1 = (l1+l2)-1;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
402
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
403 assert(NULL!=aln_param);
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
404
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
405 /*
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
406 Log(&rLog, LOG_FORCED_DEBUG, "DNAFLAG=%d seq_no=%d l1=%d l2=%d window=%d ktup=%d signif=%d wind_gap=%d",
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
407 DNAFLAG, seq_no, l1, l2, window, ktup, signif,
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
408 wind_gap);
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
409 */
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
410
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
411 slopes = (char *) CKCALLOC(tl1+1, sizeof(char));
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
412 displ = (int *) CKCALLOC(tl1+1, sizeof(int));
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
413 diag_index = (int *) CKMALLOC((tl1+1) * sizeof(int));
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
414
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
415 for (i=1; i<=tl1; ++i) {
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
416 /* unnecessary, because we calloced: slopes[i] = displ[i] = 0; */
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
417 diag_index[i] = i;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
418 }
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
419
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
420 for (i=1;i<=aln_param->ktup;i++)
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
421 pot[i] = (int) pow((double)(max_res_code+1),(double)(i-1));
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
422 limit = (int) pow((double)(max_res_code+1),(double)aln_param->ktup);
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
423
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
424
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
425
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
426 /* increment diagonal score for each k_tuple match */
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
427
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
428 for (i=1; i<=limit; ++i) {
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
429 vn1=zzc[i];
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
430 while (TRUE) {
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
431 if (!vn1) break;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
432 vn2 = zzd[i];
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
433 while (0 != vn2) {
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
434 osptr = vn1-vn2+l2;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
435 ++displ[osptr];
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
436 vn2 = zzb[vn2];
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
437 }
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
438 vn1=zza[vn1];
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
439 }
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
440 }
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
441
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
442
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
443 /* choose the top SIGNIF diagonals
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
444 */
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
445
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
446 #ifdef QSORT_REPLACEMENT
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
447 /* This was an attempt to replace des_quick_sort with qsort(),
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
448 * which turns out to be much slower, so don't use this
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
449 */
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
450
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
451 /* FIXME: if we use this branch, we don't need to init diag_index
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
452 * before, because that is done in QSortAndTrackIndex()
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
453 * automatically.
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
454 */
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
455 #if 0
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
456 for (i=1; i<=tl1; i++) {
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
457 Log(&rLog, LOG_FORCED_DEBUG, "b4 sort disp[%d]=%d diag_index[%d]=%d", i, diag_index[i], i, displ[i]);
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
458 }
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
459 #endif
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
460
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
461 QSortAndTrackIndex(&(diag_index[1]), &(displ[1]), tl1, 'a', TRUE);
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
462
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
463 #if 0
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
464 for (i=1; i<=tl1; i++) {
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
465 Log(&rLog, LOG_FORCED_DEBUG, "after sort disp[%d]=%d diag_index[%d]=%d", i, diag_index[i], i, displ[i]);
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
466 }
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
467 #endif
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
468
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
469 #else
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
470
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
471 des_quick_sort(displ, diag_index, tl1);
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
472
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
473 #endif
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
474
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
475 j = tl1 - aln_param->signif + 1;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
476
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
477 if (j < 1) {
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
478 j = 1;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
479 }
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
480
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
481 /* flag all diagonals within WINDOW of a top diagonal */
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
482
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
483 for (i=tl1; i>=j; i--) {
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
484 if (displ[i] > 0) {
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
485 pos = diag_index[i];
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
486 l = (1 > pos - aln_param->window) ?
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
487 1 : pos - aln_param->window;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
488 m = (tl1 < pos + aln_param->window) ?
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
489 tl1 : pos + aln_param->window;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
490 for (; l <= m; l++)
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
491 slopes[l] = 1;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
492 }
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
493 }
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
494
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
495 for (i=1; i<=tl1; i++) {
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
496 displ[i] = 0;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
497 }
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
498
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
499 curr_frag=*maxsf=0;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
500
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
501 for (i=1; i<=(l1-aln_param->ktup+1); ++i) {
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
502 encrypt=flag=0;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
503 for (j=1; j<=aln_param->ktup; ++j) {
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
504 residue = seq_array[seq_no][i+j-1];
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
505 if ((residue<0) || (residue>max_res_code)) {
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
506 flag=TRUE;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
507 break;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
508 }
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
509 encrypt += ((residue)*pot[j]);
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
510 }
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
511 if (flag) {
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
512 continue;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
513 }
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
514 ++encrypt;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
515
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
516 vn2=zzd[encrypt];
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
517
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
518 flag=FALSE;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
519 while (TRUE) {
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
520 if (!vn2) {
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
521 flag=TRUE;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
522 break;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
523 }
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
524 osptr=i-vn2+l2;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
525 if (1 != slopes[osptr]) {
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
526 vn2=zzb[vn2];
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
527 continue;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
528 }
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
529 flen=0;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
530 fs=aln_param->ktup;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
531 next=*maxsf;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
532
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
533 /*
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
534 * A-loop
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
535 */
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
536
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
537 while (TRUE) {
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
538 if (!next) {
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
539 ++curr_frag;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
540 if (curr_frag >= 2*max_aln_length) {
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
541 Log(&rLog, LOG_VERBOSE, "(Partial alignment)");
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
542 goto free_and_exit; /* Yesss! Always wanted to
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
543 * use a goto (AW) */
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
544 }
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
545 displ[osptr]=curr_frag;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
546 put_frag(fs, i, vn2, flen, curr_frag, &next, maxsf, accum);
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
547
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
548 } else {
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
549 tv1=accum[1][next];
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
550 tv2=accum[2][next];
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
551
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
552 if (frag_rel_pos(i, vn2, tv1, tv2, aln_param->ktup)) {
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
553 if (i-vn2 == accum[1][next]-accum[2][next]) {
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
554 if (i > accum[1][next]+(aln_param->ktup-1)) {
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
555 fs = accum[0][next]+aln_param->ktup;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
556 } else {
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
557 rmndr = i-accum[1][next];
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
558 fs = accum[0][next]+rmndr;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
559 }
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
560 flen=next;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
561 next=0;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
562 continue;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
563
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
564 } else {
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
565 if (0 == displ[osptr]) {
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
566 subt1=aln_param->ktup;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
567 } else {
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
568 if (i > accum[1][displ[osptr]]+(aln_param->ktup-1)) {
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
569 subt1=accum[0][displ[osptr]]+aln_param->ktup;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
570 } else {
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
571 rmndr=i-accum[1][displ[osptr]];
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
572 subt1=accum[0][displ[osptr]]+rmndr;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
573 }
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
574 }
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
575 subt2=accum[0][next] - aln_param->wind_gap + aln_param->ktup;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
576 if (subt2>subt1) {
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
577 flen=next;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
578 fs=subt2;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
579 } else {
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
580 flen=displ[osptr];
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
581 fs=subt1;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
582 }
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
583 next=0;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
584 continue;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
585 }
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
586 } else {
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
587 next=accum[4][next];
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
588 continue;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
589 }
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
590 }
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
591 break;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
592 }
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
593 /*
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
594 * End of Aloop
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
595 */
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
596
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
597 vn2=zzb[vn2];
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
598 }
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
599 }
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
600
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
601 free_and_exit:
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
602 CKFREE(displ);
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
603 CKFREE(slopes);
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
604 CKFREE(diag_index);
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
605
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
606 return;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
607 }
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
608 /* end of pair_align */
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
609
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
610
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
611
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
612 /**
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
613 *
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
614 * Will compute ktuple scores and store in tmat
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
615 * Following values will be set: tmat[i][j], where
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
616 * istart <= i <iend
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
617 * and
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
618 * jstart <= j < jend
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
619 * i.e. zero-offset
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
620 * tmat data members have to be preallocated
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
621 *
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
622 * if ktuple_param_t *aln_param == NULL defaults will be used
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
623 */
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
624 void
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
625 KTuplePairDist(symmatrix_t *tmat, mseq_t *mseq,
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
626 int istart, int iend,
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
627 int jstart, int jend,
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
628 ktuple_param_t *param_override,
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
629 progress_t *prProgress,
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
630 unsigned long int *ulStepNo, unsigned long int ulTotalStepNo)
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
631 {
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
632 /* this first group of variables were previously static
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
633 and hence un-parallelisable */
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
634 char **seq_array;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
635 int maxsf;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
636 int **accum;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
637 int max_aln_length;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
638 /* divide score with length of smallest sequence */
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
639 int *zza, *zzb, *zzc, *zzd;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
640 int private_step_no = 0;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
641
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
642 int i, j, dsr;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
643 double calc_score;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
644 int max_res_code = -1;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
645
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
646 int max_seq_len;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
647 int *seqlen_array;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
648 /* progress_t *prProgress; */
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
649 /* int uStepNo, uTotalStepNo; */
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
650 ktuple_param_t aln_param = default_protein_param;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
651 bool bPrintCR = (rLog.iLogLevelEnabled<=LOG_VERBOSE) ? FALSE : TRUE;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
652
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
653
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
654 if(prProgress == NULL) {
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
655 NewProgress(&prProgress, LogGetFP(&rLog, LOG_INFO),
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
656 "Ktuple-distance calculation progress", bPrintCR);
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
657 }
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
658
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
659 /* conversion to old style data types follows
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
660 *
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
661 */
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
662
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
663 seqlen_array = (int*) CKMALLOC((mseq->nseqs+1) * sizeof(int));
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
664 for (i=0; i<mseq->nseqs; i++) {
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
665 seqlen_array[i+1] = mseq->sqinfo[i].len;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
666 }
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
667
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
668 /* setup alignment parameters
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
669 */
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
670 if (SEQTYPE_PROTEIN == mseq->seqtype) {
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
671 DNAFLAG = FALSE;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
672 max_res_code = strlen(AMINO_ACID_CODES)-2;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
673 aln_param = default_protein_param;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
674
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
675 } else if (SEQTYPE_RNA == mseq->seqtype || SEQTYPE_DNA == mseq->seqtype) {
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
676 DNAFLAG = TRUE;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
677 max_res_code = strlen(NUCLEIC_ACID_CODES)-2;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
678 aln_param = default_dna_param;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
679
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
680 } else {
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
681 Log(&rLog, LOG_FATAL, "Internal error in %s: Unknown sequence type.", __FUNCTION__);
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
682 }
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
683
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
684 if (NULL!=param_override) {
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
685 aln_param.ktup = param_override->ktup;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
686 aln_param.wind_gap = param_override->wind_gap;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
687 aln_param.signif = param_override->signif;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
688 aln_param.window = param_override->window;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
689 }
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
690
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
691 /*LOG_DEBUG("DNAFLAG = %d max_res_code = %d", DNAFLAG, max_res_code);*/
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
692
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
693 /* convert mseq to clustal's old-style int encoded sequences (unit-offset)
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
694 */
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
695 max_aln_length = 0;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
696 max_seq_len = 0;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
697 seq_array = (char **) CKMALLOC((mseq->nseqs+1) * sizeof(char *));
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
698 seq_array[0] = NULL;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
699 /* FIXME check that non of the seqs is smaller than ktup (?).
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
700 * Otherwise segfault occurs
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
701 */
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
702 for (i=0; i<mseq->nseqs; i++) {
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
703 seq_array[i+1] = (char *) CKMALLOC((seqlen_array[i+1]+2) * sizeof (char));;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
704 }
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
705 for (i=0; i<mseq->nseqs; i++) {
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
706 /*LOG_DEBUG("calling encode with seq_array[%d+1] len=%d and seq=%s",
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
707 i, seqlen_array[i+1], mseq->seq[i]);*/
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
708 if (TRUE == DNAFLAG) {
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
709 encode(&(mseq->seq[i][-1]), seq_array[i+1],
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
710 seqlen_array[i+1], NUCLEIC_ACID_CODES);
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
711 } else {
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
712 encode(&(mseq->seq[i][-1]), seq_array[i+1],
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
713 seqlen_array[i+1], AMINO_ACID_CODES);
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
714 }
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
715
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
716 if (seqlen_array[i+1]>max_seq_len) {
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
717 max_seq_len = seqlen_array[i+1];
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
718 }
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
719 }
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
720 max_aln_length = max_seq_len * 2;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
721 /* see sequence.c in old source */
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
722
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
723 /* FIXME: short sequences can cause seg-fault
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
724 * because max_aln_length can get shorter
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
725 * than (max_res_code+1)^k
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
726 * FS, r222->r223 */
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
727 max_aln_length = max_aln_length > pow((max_res_code+1), aln_param.ktup)+1 ?
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
728 max_aln_length : pow((max_res_code+1), aln_param.ktup)+1;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
729
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
730 /*
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
731 *
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
732 * conversion to old style clustal done (in no time) */
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
733
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
734
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
735 accum = (int **) CKCALLOC(5, sizeof (int *));
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
736 for (i=0;i<5;i++) {
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
737 accum[i] = (int *) CKCALLOC((2*max_aln_length+1), sizeof(int));
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
738 }
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
739 zza = (int *) CKCALLOC( (max_aln_length+1), sizeof(int));
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
740 zzb = (int *) CKCALLOC( (max_aln_length+1), sizeof(int));
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
741 zzc = (int *) CKCALLOC( (max_aln_length+1), sizeof(int));
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
742 zzd = (int *) CKCALLOC( (max_aln_length+1), sizeof(int));
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
743
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
744 /* estimation of total number of steps (if istart and jstart are
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
745 * both 0) (now handled in the calling routine)
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
746 */
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
747 /* uTotalStepNo = iend*jend - iend*iend/2 + iend/2;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
748 uStepNo = 0; */
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
749 /*LOG_DEBUG("istart=%d iend=%d jstart=%d jend=%d", istart, iend, jstart, jend);*/
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
750
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
751 for (i=istart+1; i<=iend; ++i) {
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
752 /* by definition a sequence compared to itself should give
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
753 a score of 0. AW */
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
754 SymMatrixSetValue(tmat, i-1, i-1, 0.0);
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
755 make_ptrs(zza, zzc, i, seqlen_array[i], aln_param.ktup, max_res_code, seq_array);
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
756
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
757 #ifdef HAVE_OPENMP
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
758 #pragma omp critical(ktuple)
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
759 #endif
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
760 {
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
761 ProgressLog(prProgress, *ulStepNo, ulTotalStepNo, FALSE);
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
762 }
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
763
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
764 for (j=MAX(i+1, jstart+1); j<=jend; ++j) {
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
765 (*ulStepNo)++;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
766 private_step_no++;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
767 /*LOG_DEBUG("comparing pair %d:%d", i, j);*/
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
768
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
769 make_ptrs(zzb, zzd, j, seqlen_array[j], aln_param.ktup, max_res_code, seq_array);
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
770 pair_align(i, seqlen_array[i], seqlen_array[j], max_res_code, &aln_param,
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
771 seq_array, &maxsf, accum, max_aln_length, zza, zzb, zzc, zzd);
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
772
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
773 if (!maxsf) {
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
774 calc_score=0.0;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
775 } else {
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
776 calc_score=(double)accum[0][maxsf];
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
777 if (percent) {
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
778 dsr=(seqlen_array[i]<seqlen_array[j]) ?
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
779 seqlen_array[i] : seqlen_array[j];
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
780 calc_score = (calc_score/(double)dsr) * 100.0;
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
781 }
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
782 }
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
783
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
784 /* printf("%d %d %d\n", i-1, j-1, (100.0 - calc_score)/100.0); */
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
785 SymMatrixSetValue(tmat, i-1, j-1, (100.0 - calc_score)/100.0);
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
786
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
787 /* the function allows you not to compute the full matrix.
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
788 * here we explicitely make the resulting matrix a
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
789 * rectangle, i.e. we always set full rows. in other
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
790 * words, if we don't complete the full matrix then we
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
791 * don't have a full symmetry. so only use the defined
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
792 * symmetric part. AW
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
793 */
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
794 /*LOG_DEBUG("setting %d : %d = %f", j, i, tmat[i][j]);*/
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
795 /* not needed anymore since we use symmatrix_t
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
796 if (j<=iend) {
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
797 tmat[j][i] = tmat[i][j];
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
798 }
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
799 */
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
800 #ifdef HAVE_OPENMP
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
801 #pragma omp critical(ktuple)
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
802 #endif
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
803 {
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
804 Log(&rLog, LOG_DEBUG, "K-tuple distance for sequence pair %d:%d = %lg",
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
805 i, j, SymMatrixGetValue(tmat, i-1, j-1));
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
806 }
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
807 }
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
808 }
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
809 /*
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
810 Log(&rLog, LOG_FORCED_DEBUG, "uTotalStepNo=%d for istart=%d iend=%d jstart=%d jend=%d", uStepNo, istart, iend, jstart, jend);
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
811 Log(&rLog, LOG_FORCED_DEBUG, "Fabian = %d", iend*jend - iend*iend/2 + iend/2);
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
812 */
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
813
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
814 /* printf("\n\n%d\t%d\t%d\t%d\n\n", omp_get_thread_num(), uStepNo, istart, iend); */
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
815
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
816 for (i=0;i<5;i++) {
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
817 CKFREE(accum[i]);
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
818 }
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
819 CKFREE(accum);
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
820
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
821 #ifdef HAVE_OPENMP
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
822 #pragma omp critical(ktuple)
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
823 #if 0
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
824 {
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
825 printf("steps: %d\n", private_step_no);
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
826 }
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
827 #endif
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
828 #endif
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
829
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
830 CKFREE(zza);
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
831 CKFREE(zzb);
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
832 CKFREE(zzc);
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
833 CKFREE(zzd);
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
834
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
835 free(seqlen_array);
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
836
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
837 for (i=1; i<=mseq->nseqs; i++) {
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
838 CKFREE(seq_array[i]);
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
839 }
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
840 CKFREE(seq_array);
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
841 }
bc707542e5de Uploaded
clustalomega
parents:
diff changeset
842 /* end of KTuplePairDist */