« back to all changes in this revision
Viewing changes to kalign2_main.c
- browse files at revision 3
- compare with another revision
- compare with revision 3
- download diff
- download tarball
- view history from revision 3
- Committer: Bazaar Package Importer
- Author(s): Charles Plessy
- Date: 2007-02-21 09:07:33 UTC
- mfrom: (1.1.1 upstream)
- Revision ID: james.westby@ubuntu.com-20070221090733-ttdj0bl7yes4to8f
Tags: 2.03-1
* New upstream release.
* Mended the watch file.
* Updated the manpage.
* Mended the watch file.
* Updated the manpage.
- files added:
- Makefile.in
- files removed:
- Makefile
- files modified:
- README
added
removed
kalign2_main.c
1
/*
2
kalign2_main.c
3
4
Released under GPL - see the 'COPYING' file
5
6
Copyright (C) 2006 Timo Lassmann <timolassmann@gmail.com>
7
8
This program is free software; you can redistribute it and/or modify
9
it under the terms of the GNU General Public License as published by
10
the Free Software Foundation; either version 2 of the License, or
11
any later version.
12
13
This program is distributed in the hope that it will be useful,
14
but WITHOUT ANY WARRANTY; without even the implied warranty of
15
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16
GNU General Public License for more details.
17
18
You should have received a copy of the GNU General Public License
19
along with this program; if not, write to the Free Software
20
Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
21
22
Please send bug reports, comments etc. to:
23
timolassmann@gmail.com
24
*/
25
26
27
#include "kalign2.h"
28
29
30
31
unsigned int numseq = 0;
32
unsigned int numprofiles = 0;
33
float gpo = 0;
34
float gpe = 0;
35
float tgpe = 0;
36
37
int main(int argc,char **argv)
38
{
39
int i;
40
int* tree = 0;
41
int a, b, c;
42
43
struct alignment* aln = 0;
44
struct parameters* param = 0;
45
struct aln_tree_node* tree2 = 0;
46
47
param = malloc(sizeof(struct parameters));
48
49
param = interface(param,argc,argv);
50
51
aln = detect_and_read_sequences(aln,param);
52
53
if(param->ntree > numseq){
54
param->ntree = numseq;
55
}
56
57
//DETECT DNA
58
if(param->dna == -1){
59
for (i = 0; i < numseq;i++){
60
param->dna = byg_detect(aln->s[i],aln->sl[i]);
61
if(param->dna){
62
break;
63
}
64
}
65
}
66
//param->dna = 0;
67
//fprintf(stderr,"DNA:%d\n",param->dna);
68
//exit(0);
69
70
if(param->dna == 1){
71
//brief sanity check...
72
for (i = 0; i < numseq;i++){
73
if(aln->sl[i] < 6){
74
fprintf(stderr,"Dna/Rna alignments are only supported for sequences longer than 6.");
75
free(param);
76
free_aln(aln);
77
exit(0);
78
}
79
}
80
aln = make_dna(aln);
81
}
82
83
int j;
84
85
if(param->reformat){
86
for (i = 0 ;i < numseq;i++){
87
aln->nsip[i] = i;
88
for (j = 0; j < aln->sl[i];j++){
89
aln->s[i][j] = 0;
90
}
91
}
92
param->format = "fasta";//param->reformat;
93
output(aln,param);
94
exit(1);
95
}
96
97
98
99
//fast distance calculation;
100
float** submatrix = 0;
101
submatrix = read_matrix(submatrix,param); // sets gap penalties as well.....
102
103
if(!param->quiet){
104
parameter_message(param);
105
}
106
107
if(byg_start(param->alignment_type,"profPROFprofilePROFILE") != -1){
108
profile_alignment_main(aln,param,submatrix);
109
}
110
111
float** dm = 0;
112
if(param->ntree > 1){
113
if(byg_start(param->distance,"pairclustalPAIRCLUSTAL") != -1){
114
if(byg_start(param->tree,"njNJ") != -1){
115
dm = protein_pairwise_alignment_distance(aln,dm,param,submatrix,1);
116
}else{
117
dm = protein_pairwise_alignment_distance(aln,dm,param,submatrix,0);
118
}
119
}else if(byg_start("wu",param->alignment_type) != -1){
120
dm = protein_wu_distance2(aln,dm,param);
121
// param->feature_type = "wumanber";
122
}else if(param->dna == 1){
123
if(byg_start(param->tree,"njNJ") != -1){
124
dm = dna_distance(aln,dm,param,1);
125
}else{
126
dm = dna_distance(aln,dm,param,0);
127
}
128
}else{
129
if(byg_start(param->tree,"njNJ") != -1){
130
dm = protein_wu_distance(aln,dm,param,1);
131
}else{
132
dm = protein_wu_distance(aln,dm,param,0);
133
}
134
}
135
/*int j;
136
for (i = 0; i< numseq;i++){
137
for (j = 0; j< numseq;j++){
138
fprintf(stderr,"%f ",dm[i][j]);
139
}
140
fprintf(stderr,"\n");
141
}*/
142
143
if(byg_start(param->tree,"njNJ") != -1){
144
tree2 = real_nj(dm,param->ntree);
145
}else{
146
tree2 = real_upgma(dm,param->ntree);
147
}
148
if(param->print_tree){
149
print_tree(tree2,aln,param->print_tree);
150
}
151
}
152
153
tree = malloc(sizeof(int)*(numseq*3+1));
154
for ( i = 1; i < (numseq*3)+1;i++){
155
tree[i] = 0;
156
}
157
tree[0] = 1;
158
159
if(param->ntree < 2){
160
tree[0] = 0;
161
tree[1] = 1;
162
163
c = numseq;
164
tree[2] = c;
165
a = 2;
166
for ( i = 3; i < (numseq-1)*3;i+=3){
167
tree[i] = a;
168
tree[i+1] = c;
169
c++;
170
tree[i+2] = c;
171
a++;
172
}
173
}else if(param->ntree > 2){
174
ntreeify(tree2,param->ntree);
175
}else{
176
tree = readtree(tree2,tree);
177
for (i = 0; i < (numseq*3);i++){
178
tree[i] = tree[i+1];
179
}
180
free(tree2->links);
181
free(tree2->internal_lables);
182
free(tree2);
183
}
184
185
186
187
188
//get matrices...
189
struct feature_matrix* fm = 0;
190
191
struct ntree_data* ntree_data = 0;
192
193
int** map = 0;
194
if(param->ntree > 2){
195
ntree_data = malloc(sizeof(struct ntree_data));
196
ntree_data->realtree = tree2;
197
ntree_data->aln = aln;
198
ntree_data->profile = 0;
199
ntree_data->map = 0;
200
ntree_data->ntree = param->ntree;
201
ntree_data->submatrix = submatrix;
202
ntree_data->tree = tree;
203
204
ntree_data = ntree_alignment(ntree_data);
205
map = ntree_data->map;
206
tree = ntree_data->tree;
207
for (i = 0; i < (numseq*3);i++){
208
tree[i] = tree[i+1];
209
}
210
free(ntree_data);
211
}else if (param->feature_type){
212
fm = get_feature_matrix(fm,aln,param->feature_type);
213
if(!fm){
214
for (i = 32;i--;){
215
free(submatrix[i]);
216
}
217
free(submatrix);
218
free_param(param);
219
free(map);
220
free(tree);
221
exit(0);
222
}
223
224
map = feature_hirschberg_alignment(aln,tree,submatrix,map,fm);
225
//exit(0);
226
//map = feature_alignment(aln,tree,submatrix, map,fm);
227
228
}else if (byg_start("fast",param->alignment_type) != -1){
229
map = default_alignment(aln,tree,submatrix, map);
230
}else if(param->dna == 1){
231
map = dna_alignment(aln,tree,submatrix, map);
232
/*}else if (byg_start("test",param->alignment_type) != -1){
233
map = test_alignment(aln,tree,submatrix, map,param->internal_gap_weight,param->smooth_window,param->smooth_strength);
234
}else if (param->aa){
235
map = aa_alignment(aln,tree,submatrix, map,param->aa);
236
}else if (param->alter_gaps){
237
map = alter_gaps_alignment(aln,tree,submatrix,map,param->alter_gaps,param->alter_range,param->alter_weight);
238
}else if (byg_start("altergaps",param->alignment_type) != -1){
239
map = alter_gaps_alignment(aln,tree,submatrix,map,param->alter_gaps,param->alter_range,param->alter_weight);
240
}else if(byg_start("simple",param->alignment_type) != -1){
241
map = simple_hirschberg_alignment(aln,tree,submatrix, map);*/
242
}else if(byg_start("advanced",param->alignment_type) != -1){
243
map = advanced_hirschberg_alignment(aln,tree,submatrix, map,param->smooth_window,param->smooth_strength,param->internal_gap_weight);
244
}else{
245
map = hirschberg_alignment(aln,tree,submatrix, map,param->smooth_window,param->smooth_strength);
246
}
247
248
249
//clear up sequence array to be reused as gap array....
250
int *p = 0;
251
for (i = 0; i < numseq;i++){
252
p = aln->s[i];
253
for (a = 0; a < aln->sl[i];a++){
254
p[a] = 0;
255
}
256
}
257
//clear up
258
259
for (i = 0; i < (numseq-1)*3;i +=3){
260
a = tree[i];
261
b = tree[i+1];
262
aln = make_seq(aln,a,b,map[tree[i+2]]);
263
}
264
265
//for (i = 0; i < numseq;i++){
266
// fprintf(stderr,"%s %d\n",aln->sn[i],aln->nsip[i]);
267
//}
268
269
270
for (i = 0; i < numseq;i++){
271
aln->nsip[i] = 0;
272
}
273
274
275
276
aln = sort_sequences(aln,tree,param->sort);
277
278
//for (i = 0; i < numseq;i++){
279
// fprintf(stderr,"%d %d %d\n",i,aln->nsip[i],aln->sip[i][0]);
280
//}
281
282
283
output(aln,param);
284
/* if(!param->format){
285
fasta_output(aln,param->outfile);
286
}else{
287
if (byg_start("msf",param->format) != -1){
288
msf_output(aln,param->outfile);
289
}else if (byg_start("clustal",param->format) != -1){
290
clustal_output(aln,param->outfile);
291
}else if (byg_start("macsim",param->format) != -1){
292
macsim_output(aln,param->outfile,param->infile[0]);
293
}
294
}
295
free_param(param);*/
296
297
free(map);
298
free(tree);
299
return 0;
300
}
301
302
303
Loggerhead is a web-based interface for Breezy
Version: 2.0.1