1 /* -*- mode: c; tab-width: 4; c-basic-offset: 4; indent-tabs-mode: nil -*- */
3 /*********************************************************************
4 * Clustal Omega - Multiple sequence alignment
6 * Copyright (C) 2010 University College Dublin
8 * Clustal-Omega is free software; you can redistribute it and/or
9 * modify it under the terms of the GNU General Public License as
10 * published by the Free Software Foundation; either version 2 of the
11 * License, or (at your option) any later version.
13 * This file is part of Clustal-Omega.
15 ********************************************************************/
18 * RCS $Id: hhalignment-C.h 236 2011-04-14 11:30:04Z fabian $
23 * Changelog: Michael Remmert made changes to hhalign stand-alone code
24 * FS implemented some of the changes on 2010-10-28 -> MR1
26 * Note: MR seems to have changed all [aijk]++ to ++[aijk],
27 * FS did not do that on 2010-10-28
32 //////////////////////////////////////////////////////////////////////////////
34 //////////////////////////////////////////////////////////////////////////////
40 #include <iostream> // cin, cout, cerr
41 #include <fstream> // ofstream, ifstream
42 #include <stdio.h> // printf
49 #include <stdlib.h> // exit
50 #include <string> // strcmp, strstr
51 #include <math.h> // sqrt, pow
52 #include <limits.h> // INT_MIN
53 #include <float.h> // FLT_MIN
54 #include <time.h> // clock
55 #include <ctype.h> // islower, isdigit etc
56 #include "util-C.h" // imax, fmax, iround, iceil, ifloor, strint, strscn, strcut, substr, uprstr, uprchr, Basename etc.
57 #include "list.h" // list data structure
58 #include "hash.h" // hash data structure
60 #include "hhutil-C.h" // imax, fmax, iround, iceil, ifloor, strint, strscn, strcut, substr, uprstr, uprchr, Basename etc.
65 enum {KEEP_NOT = 0, KEEP_CONDITIONALLY, KEEP_ALWAYS};
67 //////////////////////////////////////////////////////////////////////////////
69 //////////////////////////////////////////////////////////////////////////////
72 //////////////////////////////////////////////////////////////////////////////
74 //////////////////////////////////////////////////////////////////////////////
75 Alignment::Alignment(int maxseq, int maxres)
78 //printf(">>>>>>>>%s:%s:%d: maxseq=%d, maxres=%d\n", __FUNCTION__, __FILE__, __LINE__, maxseq, maxres); /* (FS) */
79 longname = new(char[DESCLEN]);
80 sname = new(char*[maxseq+2]); /* MR1 */
81 seq = new(char*[maxseq+2]); /* MR1 */
83 X = new(char*[maxseq+2]); /* MR1 */
84 I = new(short unsigned int*[maxseq+2]); /* MR1 */
85 keep = new(char[maxseq+2]); /* MR1 */
86 display = new(char[maxseq+2]); /* MR1 */
87 wg = new(float[maxseq+2]); /* MR1 */
88 nseqs = new(int[maxres+2]); /* MR1 */
90 nres=NULL; // number of residues per sequence k
91 first=NULL; // first residue in sequence k
92 last=NULL; // last residue in sequence k
93 ksort=NULL; // sequence indices sorted by descending nres[k]
94 name[0]='\0'; // no name defined yet
95 longname[0]='\0'; // no name defined yet
96 fam[0]='\0'; // no name defined yet
97 file[0]='\0'; // no name defined yet
98 readCommentLine = '0'; /* MR1 */
101 //////////////////////////////////////////////////////////////////////////////
103 //////////////////////////////////////////////////////////////////////////////
104 Alignment::~Alignment()
106 delete[] longname; longname = NULL;
107 for(int k=0; k<N_in; k++)
109 delete[] sname[k]; sname[k] = NULL;
110 delete[] seq[k]; seq[k] = NULL;
111 delete[] X[k]; X[k] = NULL;
112 delete[] I[k]; I[k] = NULL;
114 delete[] sname; sname = NULL;
115 delete[] seq; seq = NULL;
116 delete[] X; X = NULL;
117 delete[] I; I = NULL;
118 delete[] l; l = NULL;
119 delete[] keep; keep = NULL;
120 delete[] display; display = NULL;
121 delete[] wg; wg = NULL;
122 delete[] nseqs; nseqs = NULL;
123 delete[] nres; nres = NULL;
124 delete[] first; first = NULL;
125 delete[] last; last = NULL;
126 delete[] ksort; ksort = NULL;
131 * @brief Reads in an alignment from file into matrix seq[k][l] as ASCII
134 Alignment::Read(FILE* inf, char infile[], char* firstline)
136 int l; // Postion in alignment incl. gaps (first=1)
137 int h; // Position in input line (first=0)
138 int k; // Index of sequence being read currently (first=0)
139 char line[LINELEN]=""; // input line
140 //char cur_seq[MAXCOL]; // Sequence currently read in
141 char *cur_seq=new(char[par.maxColCnt]);
142 char* cur_name; // Sequence currently read in
143 int linenr=0; // current line number in input file
144 char skip_sequence=0;
145 RemoveExtension(file,infile); //copy rootname (w/o path) of infile into file variable of class object
147 kss_dssp=ksa_dssp=kss_pred=kss_conf=kfirst=-1;
152 cur_seq[0]=' '; // overwrite '\0' character at beginning to be able to do strcpy(*,cur_seq)
155 // Does firstline already contain first line of file?
156 if (firstline!= NULL) strcpy(line,firstline);
158 /////////////////////////////////////////////////////////////////////////
159 // Read infile line by line
160 /* FIXME: not safe to use MAXSEQ, however, don't think we ever get here (FS) */
161 while(firstline || (fgetline(line,LINELEN,inf) && (k<MAXSEQ))) /* FIXME: FS introduced () around &&, precedence! MR1 */
165 if (line[0]=='>') //line contains sequence name
169 if (v>=1 && k>=MAXSEQ)
170 cerr<<endl<<"WARNING: maximum number "<<MAXSEQ<<" of sequences exceded in file "<<infile<<"\n";
173 cur_name=line+1; //beginning of current sequence name
174 if (k>=0) //if this is at least the second name line
176 if (strlen(cur_seq)==0)
178 cerr<<endl<<"Error: sequence "<<sname[k]<<" contains no residues."<<endl;
182 // Create space for residues and paste new sequence in
183 seq[k]=new(char[strlen(cur_seq)+2]);
184 if (!seq[k]) MemoryError("array for input sequences");
185 X[k]=new(char[strlen(cur_seq)+2]);
186 if (!X[k]) MemoryError("array for input sequences");
187 I[k]=new(short unsigned int[strlen(cur_seq)+2]);
188 if (!I[k]) MemoryError("array for input sequences");
189 strcpy(seq[k],cur_seq);
194 l=1; //position in current sequence (first=1)
196 // display[k]= 0: do not show in Q-T alignments 1: show if not filtered out later 2: show in any case (do not filter out)
197 // keep[k] = 0: do not include in profile 1: include if not filtered out later 2: include in any case (do not filter out)
198 /* {KEEP_NOT=0, KEEP_CONDITIONALLY=1, KEEP_ALWAYS=2} */
199 if (line[1]=='@') cur_name++; //skip @-character in name
200 if (!strncmp(line,">ss_dssp",8)) {
201 if (kss_dssp<0) {display[k]=2; n_display++; keep[k]=KEEP_NOT; kss_dssp=k; N_ss++;} else {skip_sequence=1; k--; continue;}
203 else if (!strncmp(line,">sa_dssp",8)) {
204 if (ksa_dssp<0) {display[k]=KEEP_ALWAYS; n_display++; keep[k]=KEEP_NOT; ksa_dssp=k; N_ss++;} else {skip_sequence=1; k--; continue;}
206 else if (!strncmp(line,">ss_pred",8)) {
207 if (kss_pred<0) {display[k]=KEEP_ALWAYS; n_display++; keep[k]=KEEP_NOT; kss_pred=k; N_ss++;} else {skip_sequence=1; k--; continue;}
209 else if (!strncmp(line,">ss_conf",8)) {
210 if (kss_conf<0) {display[k]=KEEP_ALWAYS; n_display++; keep[k]=KEEP_NOT; kss_conf=k; N_ss++;} else {skip_sequence=1; k--; continue;}
212 else if (!strncmp(line,">ss_",4) || !strncmp(line,">sa_",4)) {
213 display[k]=KEEP_ALWAYS; n_display++; keep[k]=KEEP_NOT; N_ss++;
215 else if (!strncmp(line,">aa_",4)) { // ignore sequences beginning with ">aa_"
216 skip_sequence=1; k--; continue;
218 //store first real seq
222 strwrd(word,line); // Copies first word in ptr to str
223 if (strstr(word,"_consensus"))
224 {display[k]=2; keep[k]=0; n_display++; kfirst=k;} /* MR1 */
226 {display[k]=keep[k]=KEEP_ALWAYS; n_display++; kfirst=k;}
228 //store all sequences
229 else if (par.mark==0) {display[k]=keep[k]=KEEP_CONDITIONALLY; n_display++;}
230 //store sequences up to nseqdis
231 else if (line[1]=='@'&& n_display-N_ss<par.nseqdis) {display[k]=keep[k]=KEEP_ALWAYS; n_display++;}
232 else {display[k]=KEEP_NOT; keep[k]=KEEP_CONDITIONALLY;}
234 // store sequence name
235 if (v>=4) printf("Reading seq %-16.16s k=%3i n_displ=%3i display[k]=%i keep[k]=%i\n",cur_name,k,n_display,display[k],keep[k]);
236 sname[k] = new(char[strlen(cur_name)+1]);
237 if (!sname[k]) {MemoryError("array for sequence names");}
238 strcpy(sname[k],cur_name);
239 } // end if(line contains sequence name)
241 else if (line[0]=='#') // Commentary line?
243 // #PF01367.9 5_3_exonuc: 5'-3' exonuclease, C-terminal SAM fold; PDB 1taq, 1bgx (T:271-174), 1taq (271-174)
244 if (name[0]) continue; // if already name defined: skip commentary line
246 ptr1=strscn_(line+1); // set ptr1 to first non-whitespace character after '#' -> AC number
247 strncpy(longname,ptr1,DESCLEN-1); // copy whole commentary line after '# ' into longname
248 longname[DESCLEN-1]='\0';
249 strtr(longname,""," ");
250 ptr2=strcut_(ptr1); // cut after AC number and set ptr2 to first non-whitespace character after AC number
251 // strcpy(fam,ptr1); // copy AC number to fam
252 // if (!strncmp(fam,"PF",2)) strcut_(fam,'.'); // if PFAM identifier contains '.' cut it off
253 // strcut_(ptr2); // cut after first word ...
254 strcpy(name,ptr1); // ... and copy first word into name
255 readCommentLine = '1'; /* MR1 */
258 //line contains sequence residues or SS information and does not belong to a >aa_ sequence
259 else if (!skip_sequence)
261 if (v>=4) cout<<line<<"\n"; //DEBUG
264 cerr<<endl<<"WARNING: No sequence name preceding following line in "<<infile<<":\n\'"<<line<<"\'\n";
268 h=0; //counts characters in current line
270 // Check whether all characters are correct; store into cur_seq
271 if (keep[k] || (k == kfirst) ) // normal line containing residues /* MR1 */
273 while (h<LINELEN && line[h]>'\0' && l</*MAXCOL*/par.maxColCnt-1)
275 if (aa2i(line[h])>=0) // ignore white-space characters ' ', \t and \n (aa2i()==-1)
276 {cur_seq[l]=line[h]; l++;}
277 else if (aa2i(line[h])==-2 && v)
278 cerr<<endl<<"WARNING: invalid symbol \'"<<line[h]<<"\' at pos. "<<h<<" in line "<<linenr<<" of "<<infile<<"\n";
282 else if (k==kss_dssp) // lines with dssp secondary structure states (. - H E C S T G B)
284 while (h<LINELEN && line[h]>'\0' && l</*MAXCOL*/par.maxColCnt-1)
286 if (ss2i(line[h])>=0 && ss2i(line[h])<=7)
287 {cur_seq[l]=ss2ss(line[h]); l++;}
289 cerr<<endl<<"WARNING: invalid symbol \'"<<line[h]<<"\' at pos. "<<h<<" in line "<<linenr<<" of "<<infile<<"\n";
293 else if (k==ksa_dssp) // lines with dssp solvent accessibility states (. - ???)
295 while (h<LINELEN && line[h]>'\0' && l</*MAXCOL*/par.maxColCnt-1)
297 if (sa2i(line[h])>=0)
298 cur_seq[l++]=line[h];
300 cerr<<endl<<"WARNING: invalid symbol \'"<<line[h]<<"\' at pos. "<<h<<" in line "<<linenr<<" of "<<infile<<"\n";
304 else if (k==kss_pred) // lines with predicted secondary structure (. - H E C)
306 while (h<LINELEN && line[h]>'\0' && l</*MAXCOL*/par.maxColCnt-1)
308 if (ss2i(line[h])>=0 && ss2i(line[h])<=3)
309 {cur_seq[l]=ss2ss(line[h]); l++;}
311 cerr<<endl<<"WARNING: invalid symbol \'"<<line[h]<<"\' at pos. "<<h<<" in line "<<linenr<<" of "<<infile<<"\n";
315 else if (k==kss_conf) // lines with confidence values should contain only 0-9, '-', or '.'
317 while (h<LINELEN && line[h]>'\0' && l</*MAXCOL*/par.maxColCnt-1)
319 if (line[h]=='-' || line[h]=='.' || (line[h]>='0' && line[h]<='9'))
320 {cur_seq[l]=line[h]; l++;}
322 cerr<<endl<<"WARNING: invalid symbol \'"<<line[h]<<"\' at pos. "<<l<<" in line "<<linenr<<" of "<<infile<<"\n";
326 else if (display[k]) // other lines such as >sa_pred etc
328 while (h<LINELEN && line[h]>'\0' && l</*MAXCOL*/par.maxColCnt-1)
330 if (line[h]=='-' || line[h]=='.' || (line[h]>='0' && line[h]<='9') || (line[h]>='A' && line[h]<='B'))
331 {cur_seq[l]=line[h]; l++;}
333 cerr<<endl<<"WARNING: invalid symbol \'"<<line[h]<<"\' at pos. "<<l<<" in line "<<linenr<<" of "<<infile<<"\n";
337 if (v && l>=/*MAXCOL*/par.maxColCnt-1)
339 cerr<<endl<<"WARNING: maximum number of residues "<</*MAXCOL*/par.maxColCnt-2<<" exceded in sequence "<<sname[k]<<"\n";
342 cur_seq[l]='\0'; //Ensure that cur_seq ends with a '\0' character
346 /////////////////////////////////////////////////////////////////////////
349 if (k>=0) //if at least one sequence was read in
351 seq[k]=new(char[strlen(cur_seq)+2]);
352 if (!seq[k]) MemoryError("array for input sequences");
353 X[k]=new(char[strlen(cur_seq)+2]);
354 if (!X[k]) MemoryError("array for input sequences");
355 I[k]=new(short unsigned int[strlen(cur_seq)+2]);
356 if (!I[k]) MemoryError("array for input sequences");
357 strcpy(seq[k],cur_seq);
360 {cerr<<endl<<"Error: no sequences found in file "<<infile<<"\n"; exit(1);}
364 // Set name, longname, fam
365 if (!*name) // longname, name and family were not set by '#...' line yet -> extract from first sequence
368 // strtr(sname[kfirst],"~"," "); // 'transpose': replaces the tilde with a blanc everywhere in sname[kfirst]
369 strncpy(longname,sname[kfirst],DESCLEN-1); // longname is name of first sequence
370 longname[DESCLEN-1]='\0';
371 strncpy(name,sname[kfirst],NAMELEN-1); // Shortname is first word of longname...
372 name[NAMELEN-1]='\0';
373 ptr = strcut(name); // ...until first white-space character
374 if (ptr && islower(ptr[0]) && ptr[1]=='.' && isdigit(ptr[2])) //Scop family code present as second word?
376 lwrstr(name); // Transform upper case to lower case
377 strcut(ptr); // Non-white-space characters until next white-space character..
378 strcpy(fam,ptr); // ...are the SCOP familiy code
380 else if (name[0]=='P' && name[1]=='F' && isdigit(name[2]) && isdigit(name[3]) ) //Pfam code
382 strcpy(fam,name); // set family name = Pfam code
388 delete[] cur_seq; cur_seq = NULL;
390 // Checking for warning messages
392 if (v>=2) cout<<"Read "<<infile<<" with "<<N_in<<" sequences\n";
393 if (v>=3) cout<<"Query sequence for alignment has number "<<kfirst<<" (0 is first)\n";
398 * At this point GetSeqsFromHMM() slots in, however,
399 * only needed in hhbliys.C, so will skip it for moment, MR1
403 /////////////////////////////////////////////////////////////////////////////
405 * @brief Convert ASCII in seq[k][l] to int (0-20) in X[k][i],
406 * throw out all insert states, record their number in I[k][i]
407 * and store sequences to be displayed in seq[k] */
408 /////////////////////////////////////////////////////////////////////////////
410 Alignment::Compress(const char infile[])
412 int i; // Index for match state (first=1)
413 int l; // Postion in alignment incl. gaps (first=1)
414 int k; // Index for sequences (first=0)
415 int a; // amino acid index
417 int unequal_lengths=0; /* k: seq k doesn't have same number
418 of match states as seq 0 => WARNING */
419 /* points to next character in seq[k] to be written */
420 /*static short unsigned int h[MAXSEQ];*/
421 /*short*/ unsigned int *h = NULL; /* short may lead to overflow for long alignments, FS, r235 -> r236 */
423 h = new(/*short*/ unsigned int[N_in+2]); /* short -> overflow, FS, r235 -> r236 */
424 float *percent_gaps = NULL; /* FS, 2010-Nov */
425 char *match_state = NULL; /* FS, 2010-Nov */
429 for (k=0;k<N_in; k++)
435 cout<<"Using match state assignment by capital letters (a2m format)\n";
436 else if (par.M==2) cout<<"Using percentage-rule match state assignment\n";
437 else if (par.M==3) cout<<"Using residues of first sequence as match states\n";
440 // Create matrices X and I with amino acids represented by integer numbers
444 /////////////////////////////////////////////////////////////////////////
445 /* a2m/a3m format: match states capital case,
446 inserts lower case, delete states '-', inserted gaps '.'
447 The confidence values for ss prediction are interpreted as follows:
448 0-9:match states(!) '-' :match state '.':insert */
452 // Warn if alignment is ment to be -M first or -M NN instead of A2M/A3M
453 if (v>=2 && strchr(seq[kfirst],'-') ) // Seed/query sequence contains a gap ...
455 for (k=1; k<N_in; k++)
456 if (strpbrk(seq[k],"abcdefghiklmnpqrstuvwxyz.")) break;
457 if (k==N_in) // ... but alignment contains no lower case residue
458 printf("WARNING: input alignment %s looks like aligned FASTA instead of A2M/A3M format. Consider using '-M first' or '-M 50'\n",infile);
461 // Remove '.' characters from seq[k]
462 for(k=0; k<N_in; k++)
464 char* ptrS=seq[k]; // pointer to source: character in seq[k]
465 char* ptrD=seq[k]; // pointer to destination: seq[k]
466 while(1) // omit '.' symbols
468 if (*ptrS!='.') {*ptrD=*ptrS; ptrD++;} //leave out '.' symbols
473 L=/*MAXRES*/par.maxResLen-2; // needed because L=imin(L,i)
474 for (k=0; k<N_in; k++)
476 i=1; l=1; // start at i=1, not i=0!
477 if (keep[k]) //skip >ss_dssp, >ss_pred, >ss_conf, >aa_... sequences
479 while((c=seq[k][l++])) // assign residue to c at same time
481 if (c>='a' && c<='z') I[k][i-1]++;//insert state = lower case character
482 else if (c!='.') //match state = upper case character
490 else if (k==kss_dssp || k==kss_pred) // does alignment contain sequence of secondary structure states?
492 while((c=seq[k][l++])) // assign residue to c at same time
493 if (c!='.' && !(c>='a' && c<='z')) X[k][i++]=ss2i(c); //match state = upper case character
495 else if (k==ksa_dssp) // does alignment contain sequence of prediction confidence values?
497 while((c=seq[k][l++])) // assign residue to c at same time
498 if (c!='.' && !(c>='a' && c<='z')) X[k][i++]=sa2i(c); //match state = upper case character
500 else if (k==kss_conf) // does alignment contain sequence of prediction confidence values?
502 while((c=seq[k][l++])) // assign residue to c at same time
503 if (c!='.') X[k][i++]=cf2i(c); //match state = 0-9 or '-'
505 else if (k==kfirst) // does alignment contain sequence of prediction confidence values?
507 while((c=seq[k][l++])) // assign residue to c at same time
517 if (L!=i && L!=/*MAXRES*/par.maxResLen-2 && !unequal_lengths) unequal_lengths=k; //sequences have different lengths
520 if (unequal_lengths) break;
522 //Replace GAP with ENDGAP for all end gaps /* MR1 */
523 for (k=0; k<N_in; ++k)
525 if (!keep[k]) continue;
526 for (i=1; i<=L && X[k][i]==GAP; i++) X[k][i]=ENDGAP; /* MR1: NOTE i++ <- ++i */
527 for (i=L; i>=1 && X[k][i]==GAP; i--) X[k][i]=ENDGAP; /* MR1 */
530 for (i=1; i<=L; i++) this->l[i]=i; //assign column indices to match states
533 cout<<"\nError: Alignment in "<<infile<<" contains no match states. Consider using -M first or -M <int> option"<<endl;
537 if (L==/*MAXRES*/par.maxResLen-2 && v>=2)
539 printf("WARNING: Number of match columns too large. Only first %i match columns will be kept!\n",L);
542 if (v>=2) cout<<"Alignment in "<<infile<<" contains "<<L<<" match states\n";
545 /////////////////////////////////////////////////////////////////////////
546 // gap-rule assignment of match states
548 int nl[NAA+2]; //nl[a] = number of seq's with amino acid a at position l
549 /* Note: allocating statically is fine most of the time
550 but when the sequences/profiles get really long
551 we might run out of memory, so must really do it dynamically.
552 had to move declaration of float *percent_gaps out of switch()
554 //float percent_gaps[MAXCOL]; //percentage of gaps in column k (with weighted sequences)
555 percent_gaps = new(float[par.maxColCnt]);
557 //determine number of columns L in alignment
558 L=strlen(seq[kfirst])-1;
560 // Conversion to integer representation, checking for unequal lengths and initialization
561 if (nres==NULL) nres=new(int[N_in]);
562 for (k=0; k<N_in; k++)
564 if (!keep[k]) continue;
569 X[k][l]=aa2i(seq[k][l]);
570 if (X[k][l]<NAA) nr++;
573 if (seq[k][L+1]!='\0' && !unequal_lengths) unequal_lengths=k;
575 if (unequal_lengths) break;
577 // Quick and dirty calculation of the weight per sequence wg[k]
578 for (l=1; l<=L; l++) // for all positions l in alignment
580 int naa=0; //number of different amino acids
581 for (a=0; a<20; a++) nl[a]=0;
582 for (k=0; k<N_in; k++) if (keep[k]) nl[ (int)X[k][l]]++;
583 for (a=0; a<20; a++) if(nl[a]) naa++;
584 if (!naa) naa=1; //naa=0 when column consists of only gaps and Xs (=ANY)
585 for (k=0; k<N_in; k++)
586 if (keep[k] && (X[k][l]<20) )
588 //wg[k]+=1.0/float(nl[ (int)X[k][l]]*naa*nres[k]+30.0); /* original version */
589 wg[k] += 1.0/float(nl[ (int)X[k][l]]*naa*(nres[k]+30.0)); /* MR1 */
590 // wg[k] += 1.0/float(nl[ (int)X[k][l]]*(nres[k]+30.0)); /* MR1 commented out */
591 // wg[k] += (naa-1.0)/float(nl[ (int)X[k][l]]*(nres[k]+30.0)); /* MR1 commented out */
595 //Replace GAP with ENDGAP for all end gaps
596 for (k=0; k<N_in; ++k)
598 if (!keep[k]) continue;
599 for (i=1; i<=L && X[k][i]==GAP; i++) X[k][i]=ENDGAP; /* MR1: NOTE i++ <- ++i */
600 for (i=L; i>=1 && X[k][i]==GAP; i--) X[k][i]=ENDGAP; /* MR1 */
603 // Add up percentage of gaps
608 for (k=0; k< N_in; k++){
610 if ( X[k][l]<GAP) res+=wg[k]; /* MR1, AA or ANY, changed from <ANY */
611 else if ( X[k][l] != ENDGAP) gap+=wg[k]; /* MR1, else: GAP. ENDGAPs are ignored for counting percentage */
614 percent_gaps[l]=100.*gap/(res+gap);
615 if (v>=4) cout<<"percent gaps["<<l<<"]="<<percent_gaps[l]<<" first seq:"<<seq[0][l]<<"\n";
618 /* Insert states 'bloat' the HMM,
619 throwing them out 'slims' down the HMM.
620 A slimmer HMM takes less time to construct.
621 However, the marriage of Clustal and Hhalign
622 is particularly sensitive to residues
623 at the very end of the profile; these I call
624 'telomeres'. Telomeres must not be shed when
625 throwing out insert states, for the telomeres
626 we set the match threshold to 100%.
629 #define TELOMERE_LOGIC 1
630 #define TELOMERE_DYNAMIC 0
632 #define ALWAYS_ACCEPT 101.0 /* do NOT change this parameter, must be >=100,
633 make slightly bigger than 100% -- to be sure to be sure */
634 #define DEFAULT_MGAPS 100.0 /* Soeding's default is 50, omega default prior to telomere logic was 100
635 FIXME: this used to be par.Mgaps,
636 in a later version re-introduce par.Mgaps to keep this value flexible */
637 #define TELOMER_LENGTH 10 /* this parameter must be > 0 (unless DEFAULT_MGAPS=100),
638 if it is too big (L/2) then telomere logic has no effect,
639 don't think it should be changed (much) */
640 #define TELOMER_FRACTION 0.10
641 //#define HMM_MIN_LENGTH 0.923
642 #define HMM_MIN_LENGTH 0.950
643 #define FORTRAN_OFFSET 1
644 double dDefaultMgaps;
645 dDefaultMgaps = DEFAULT_MGAPS;
647 #if TELOMERE_LOGIC /* turn telomere logic on (1) or off (0) */
650 #if TELOMERE_DYNAMIC /* keep telomere length 'dynamic' */
651 iTelomereLength = TELOMER_LENGTH > (int)(L*TELOMER_FRACTION) ? TELOMER_LENGTH : (int)(L*TELOMER_FRACTION);
653 iTelomereLength = TELOMER_LENGTH;
654 #endif /* this was dynamic telomere */
655 #endif /* this was telomere logic */
657 /* if HMMs get too small (much smaller than profile length L)
658 then one is liable to get a back-tracking error.
659 So we should ensure that the DEFAULT_MGAPS parameter does NOT
660 shrink the HMM too much.
661 take percentage-gap vector, sort it, and fix dDefaultMgaps,
662 such that at least (HMM_MIN_LENGTH)*(L) are left
664 #if MGAP_LOGIC /* try to adapt Mgaps to size of final HMM */
666 float *pfPercentGaps = NULL;
667 if (NULL == (pfPercentGaps = (float *)malloc((L+1)*sizeof(float)))){
668 printf("%s:%s:%d: could not malloc %d float for sorted percent-gaps\n",
669 __FUNCTION__, __FILE__, __LINE__, L+1);
670 dDefaultMgaps = DEFAULT_MGAPS;
673 for (l = 0; l < L; l++) {
674 pfPercentGaps[l] = percent_gaps[l+FORTRAN_OFFSET];
676 qsort(pfPercentGaps, L, sizeof(float), CompFltAsc);
678 dDefaultMgaps = pfPercentGaps[(int)(HMM_MIN_LENGTH*L)];
679 if (dDefaultMgaps < DEFAULT_MGAPS){
680 //printf("Mgaps = %f <- %f\n", DEFAULT_MGAPS, dDefaultMgaps);
681 dDefaultMgaps = DEFAULT_MGAPS;
684 //printf("Mgaps = %f\n", dDefaultMgaps);
687 free(pfPercentGaps); pfPercentGaps = NULL;
690 #endif /* tried to adapt Mgaps to size of final HMM */
692 // Throw out insert states and keep only match states
694 for (k=0; k<N_in; k++) {h[k]=1; seq[k][0]='-';}
698 float fMgaps = ALWAYS_ACCEPT;
699 if ( (l < iTelomereLength) || (L-l < iTelomereLength) ){
700 /* residue is in telomere, always retain this position */
701 fMgaps = ALWAYS_ACCEPT;
704 /* FIXME: would like to put a transition phase in here,
705 where the Mgap value gradually goes down from 100 to DEFAULT_MGAPS,
706 however, may not be necessary and will make code more clunky */
709 /* position is in centre of sequence,
710 retain position if less than DEFAULT_MGAPS% gaps at this position,
711 for example, if DEFAULT_MGAPS=30 throw out if more than 30% gap.
712 conversely, if DEFAULT_MGAPS=100 throw out if more than 100% gaps,
713 which can never happen, so always retain */
714 fMgaps = dDefaultMgaps;
716 if (percent_gaps[l] <= fMgaps)
717 #else /* this was telomere logic */
718 if (percent_gaps[l]<=float(par.Mgaps))
719 #endif /* this was Soeding default */
721 if (i>=/*MAXRES*/par.maxResLen-2) {
723 printf("WARNING: Number of match columns too large. Only first %i match columns will be kept!\n",i);
728 for (k=0; k<N_in; k++)
732 seq[k][h[k]++]=MatchChr(seq[k][l]);
736 else if (k==kss_dssp || k==kss_pred)
738 seq[k][h[k]++]=MatchChr(seq[k][l]);
739 X[k][i]=ss2i(seq[k][l]);
741 else if (k==ksa_dssp)
743 seq[k][h[k]++]=MatchChr(seq[k][l]);
744 X[k][i]=sa2i(seq[k][l]);
746 else if (k==kss_conf)
748 seq[k][h[k]++]=seq[k][l];
749 X[k][i]=cf2i(seq[k][l]);
755 for (k=0; k<N_in; k++)
756 if (keep[k] && X[k][l]<GAP)
759 seq[k][h[k]++]=InsertChr(seq[k][l]);
763 for (k=0; k<N_in; k++) seq[k][h[k]]='\0';
765 //printf("%d\t%d\t%d\tN/L/M\n", N_in, L, i); /* -------- FIXME */
767 if (v>=2) cout<<"Alignment in "<<infile<<" contains "<<L<<" columns and "<<i<<" match states\n";
768 L = i; //Number of match states
770 delete[] percent_gaps; percent_gaps = NULL;
774 ////////////////////////////////////////////////////////////////////////
775 // Using residues of first sequence as match states
777 /* Note: allocating statically is fine most of the time
778 but when the sequences/profiles get really long
779 we might run out of memory, so must really do it dynamically.
780 had to move declaration of float *percent_gaps out of switch()
782 //char match_state[MAXCOL]; //1: column assigned to match state 0: insert state
783 match_state = new(char[par.maxColCnt]);
785 // Determine number of columns L in alignment
787 if (v>=3) printf("Length of first seq = %i\n",L);
788 // Check for sequences with unequal lengths
789 for (k=1; k<N_in; k++)
790 if (int(strlen(seq[k]+1))!=L) {unequal_lengths=k; break;}
791 if (unequal_lengths) break;
793 // Determine match states: seq kfirst has residue at pos l -> match state
795 if (isalpha(seq[kfirst][l])) match_state[l]=1; else match_state[l]=0;
796 // Throw out insert states and keep only match states
797 for (k=0; k<N_in; k++) {h[k]=1; seq[k][0]='-';}
801 if (match_state[l]) // does sequence 0 have residue at position l?
803 if (i>=/*MAXRES*/par.maxResLen-2) {
805 printf("WARNING: Number of match columns too large. Only first %i match columns will be kept!\n",i);
810 for (k=0; k<N_in; k++)
814 seq[k][h[k]++]=MatchChr(seq[k][l]);
815 X[k][i]=aa2i(seq[k][l]);
818 else if (k==kss_dssp || k==kss_pred)
820 seq[k][h[k]++]=MatchChr(seq[k][l]);
821 X[k][i]=ss2i(seq[k][l]);
823 else if (k==ksa_dssp)
825 seq[k][h[k]++]=MatchChr(seq[k][l]);
826 X[k][i]=sa2i(seq[k][l]);
828 else if (k==kss_conf)
830 seq[k][h[k]++]=seq[k][l];
831 X[k][i]=cf2i(seq[k][l]);
837 for (k=0; k<N_in; k++)
838 if (keep[k] && aa2i(seq[k][l])<GAP)
841 seq[k][h[k]++]=InsertChr(seq[k][l]);
845 for (k=0; k<N_in; k++) seq[k][h[k]]='\0';
847 //Replace GAP with ENDGAP for all end gaps /* MR1 */
848 for (k=0; k<N_in; ++k)
850 if (!keep[k]) continue;
851 for (i=1; i<=L && X[k][i]==GAP; i++) X[k][i]=ENDGAP; /* MR1, note i++ <- ++i */
852 for (i=L; i>=1 && X[k][i]==GAP; i--) X[k][i]=ENDGAP; /* MR1 */
855 if (v>=2) cout<<"Alignment in "<<infile<<" contains "<<L<<" columns and "<<i<<" match states\n";
856 L = i; //Number of match states
858 delete[] match_state; match_state = NULL;
862 ///////////////////////////////////////////////////////////////////////////
868 strcut(sname[unequal_lengths]);
869 cerr<<endl<<"Error: sequences in "<<infile<<" do not all have the same number of columns, \ne.g. first sequence and sequence "<<sname[unequal_lengths]<<".\n";
870 if(par.M==1) cerr<<".\nCheck input format for '-M a2m' option and consider using '-M first' or '-M 50'\n";
874 // Avert user about -cons option?
875 if (v>=2 && !par.cons)
878 if (X[kfirst][i]==GAP)
880 printf("NOTE: Use the '-cons' option to calculate a consensus sequence as first sequence of the alignment.\n");
885 //Replace GAP with ENDGAP for all end gaps
886 for (k=0; k<N_in; k++)
888 if (!keep[k]) continue;
889 for (i=1; i<=L && X[k][i]==GAP; i++) X[k][i]=ENDGAP;
890 for (i=L; i>=1 && X[k][i]==GAP; i--) X[k][i]=ENDGAP;
895 for (k=0; k<N_in; k++)
897 if (!display[k]) continue;
898 cout<<">"<<sname[k]<<"\n";
899 if (k==kss_dssp || k==kss_pred) {for (i=1; i<=L; i++) cout<<char(i2ss(X[k][i]));}
900 else if (k==kss_conf) {for (i=1; i<=L; i++) cout<<char(i2cf(X[k][i]));}
901 else if (k==ksa_dssp) {for (i=1; i<=L; i++) cout<<char(i2sa(X[k][i]));}
904 for (i=1; i<=L; i++) cout<<char(i2aa(X[k][i]));
907 if (I[k][i]==0) cout<<"-"; else if (I[k][i]>9) cout<<"X"; else cout<<I[k][i];
912 delete[](h); h = NULL;
917 * @brief Remove sequences with seq. identity larger than seqid percent
918 *(remove the shorter of two) or coverage<cov_thr
920 * FIXME: originally max_seqid is a variable that is the cutoff
921 * above which sequences are thrown out. We want to throw out sequences
922 * when building the HMM but not for display, there we want to keep all.
923 * This should be really easy, but there is some hidden stuff going on
924 * in this function, so I did a minimal-invasive change and just stuck
925 * (effectively) a hard-wired 100 instead of the variable.
926 * At a later stage we should get rid of this function alltogether
927 * as it does gobble up some time (and is quadratic in noof sequences, I think)
930 ////////////////////////////////////////////////////////////////////////////
934 Alignment::FilterForDisplay(int max_seqid, int coverage, int qid, float qsc, int N)
938 * by just returning n_display and not doing anything
939 * I think we display everything and not do any work for it
941 return n_display; /* FS, 2010-10-04*/
944 if (par.mark) return n_display;
945 char *dummy = new(char[N_in+1]);
949 if (kss_dssp>=0) display[kss_dssp]=KEEP_NOT;
950 if (ksa_dssp>=0) display[ksa_dssp]=KEEP_NOT;
951 if (kss_pred>=0) display[kss_pred]=KEEP_NOT;
952 if (kss_conf>=0) display[kss_conf]=KEEP_NOT;
953 for (seqid=imin(10,max_seqid); n_display<N && seqid<=max_seqid; seqid++)
955 for (int k=0; k<N_in; k++) dummy[k]=display[k];
956 n_display = Filter2(dummy,coverage,qid,qsc,20,seqid,0);
957 // printf("Seqid=%3i n_display=%4i\n",seqid,n_display);
961 for (int k=0; k<N_in; k++) dummy[k]=display[k];
962 n_display = Filter2(dummy,coverage,qid,qsc,20,--(--seqid),0);
965 for (int k=0; k<N_in; k++) display[k]=dummy[k];
966 if (kss_dssp>=0) {display[kss_dssp]=KEEP_CONDITIONALLY; n_display++;}
967 if (ksa_dssp>=0) {display[ksa_dssp]=KEEP_CONDITIONALLY; n_display++;}
968 if (kss_pred>=0) {display[kss_pred]=KEEP_CONDITIONALLY; n_display++;}
969 if (kss_conf>=0) {display[kss_conf]=KEEP_CONDITIONALLY; n_display++;}
970 delete[] dummy; dummy = NULL;
974 /////////////////////////////////////////////////////////////////////////////////////
975 // Remove sequences with seq. identity larger than seqid percent (remove the shorter of two) or coverage<cov_thr
976 /////////////////////////////////////////////////////////////////////////////////////
977 inline int Alignment::Filter(int max_seqid, int coverage, int qid, float qsc, int N)
979 return Filter2(keep,coverage,qid,qsc,20,max_seqid,N);
982 /////////////////////////////////////////////////////////////////////////////
984 * @brief Select set of representative sequences in the multiple sequence alignment
987 * Remove sequences with coverage of query less than "coverage" percent
988 * Remove sequences with sequence identity to query of less than "qid" percent
989 * If Ndiff==0, remove sequences with seq. identity larger than seqid2(=max_seqid) percent
990 * If Ndiff>0, remove sequences with minimum-sequence-identity filter of between seqid1
991 * and seqid2 (%), where the minimum seqid threshold is determined such that,
992 * in all column blocks of at least WMIN=25 residues, at least Ndiff sequences are left.
993 * This ensures that in multi-domain proteins sequences covering one domain are not
994 * removed completely because sequences covering other domains are more diverse.
996 * Allways the shorter of two compared sequences is removed (=> sort sequences by length first).
997 * Please note: sequence identity of sequence x with y when filtering x is calculated as
998 * number of residues in sequence x that are identical to an aligned residue in y / number of residues in x
999 * Example: two sequences x and y are 100% identical in their overlapping region but one overlaps by 10% of its
1000 * length on the left and the other by 20% on the right. Then x has 10% seq.id with y and y has 20% seq.id. with x.
1002 //////////////////////////////////////////////////////////////////////////////
1004 Alignment::Filter2(char keep[], int coverage, int qid, float qsc, int seqid1, int seqid2, int Ndiff)
1006 // In the beginnning, keep[k] is 1 for all regular amino acid sequences and 0 for all others (ss_conf, ss_pred,...)
1007 // In the end, keep[k] will be 1 for all regular representative sequences kept in the alignment, 0 for all others
1008 char* in=new(char[N_in+1]); // in[k]=1: seq k has been accepted; in[k]=0: seq k has not yet been accepted at current seqid
1009 char* inkk=new(char[N_in+1]); // inkk[k]=1 iff in[ksort[k]]=1 else 0;
1010 int* Nmax=new(int[L+2]); // position-dependent maximum-sequence-identity threshold for filtering /* MR1, used to be called idmax*/
1011 int* idmaxwin=new(int[L+2]); // minimum value of idmax[i-WFIL,i+WFIL]
1012 int* seqid_prev=new(int[N_in+1]); // maximum-sequence-identity threshold used in previous round of filtering (with lower seqid)
1013 int* N=new(int[L+2]); // N[i] number of already accepted sequences at position i
1014 const int WFIL=25; // see previous line
1016 int diffNmax=Ndiff; // current maximum difference of Nmax[i] and Ndiff /* MR1 */
1017 int diffNmax_prev=0; // previous maximum difference of Nmax[i] and Ndiff /* MR1 */
1019 int seqid; // current maximum value for the position-dependent maximum-sequence-identity thresholds in idmax[]
1020 int seqid_step=0; // previous increment of seqid /* MR1 */
1022 float diff_min_frac; // minimum fraction of differing positions between sequence j and k needed to accept sequence k
1023 float qdiff_max_frac=0.9999-0.01*qid; // maximum allowable number of residues different from query sequence
1024 int diff; // number of differing positions between sequences j and k (counted so far)
1025 int diff_suff; // number of differing positions between sequences j and k that would be sufficient
1026 int qdiff_max; // maximum number of residues required to be different from query
1027 int cov_kj; // upper limit of number of positions where both sequence k and j have a residue
1028 int first_kj; // first non-gap position in sequence j AND k
1029 int last_kj; // last non-gap position in sequence j AND k
1030 int kk, jj; // indices for sequence from 1 to N_in
1031 int k, j; // kk=ksort[k], jj=ksort[j]
1032 int i; // counts residues
1033 int n; // number of sequences accepted so far
1037 for (n=k=0; k<N_in; k++) if (keep[k]==KEEP_ALWAYS) {in[k]=2/*KEEP_ALWAYS??*/; n++;} else in[k]=0;
1039 // Determine first[k], last[k]?
1042 first=new(int[N_in]);// first non-gap position in sequence k
1043 last =new(int[N_in]);// last non-gap position in sequence k
1044 for (k=0; k<N_in; k++) // do this for ALL sequences, not only those with in[k]==1 (since in[k] may be display[k])
1046 for (i=1; i<=L; i++) if (X[k][i]<NAA) break;
1048 for (i=L; i>=1; i--) if (X[k][i]<NAA) break;
1053 // Determine number of residues nres[k]?
1054 if ( (nres==NULL) || (sizeof(nres)<N_in*sizeof(int)) )
1056 nres=new(int[N_in]);
1057 for (k=0; k<N_in; k++) // do this for ALL sequences, not only those with in[k]==1 (since in[k] may be display[k])
1060 for (i=first[k]; i<=last[k]; i++)
1061 if (X[k][i]<NAA) nr++;
1063 // printf("%20.20s nres=%3i first=%3i last=%3i\n",sname[k],nr,first[k],last[k]);
1067 // Sort sequences according to length; afterwards, nres[ksort[kk]] is sorted by size
1070 ksort=new(int[N_in]); // never reuse alignment object for new alignment with more sequences
1071 for (k=0; k<N_in; k++) ksort[k]=k;
1072 QSortInt(nres,ksort,kfirst+1,N_in-1,-1); //Sort sequences after kfirst (query) in descending order
1074 for (kk=0; kk<N_in; kk++) inkk[kk]=in[ksort[kk]];
1076 // Initialize N[i], idmax[i], idprev[i]
1077 for (i=1; i<first[kfirst]; i++) N[i]=0;
1078 for (i=first[kfirst]; i<=last[kfirst]; i++) N[i]=1;
1079 for (i=last[kfirst]+1; i<=L; i++) N[i]=0;
1080 //for (i=1; i<=L; i++) {idmax[i]=seqid1; idmaxwin[i]=-1;}
1081 for (i=1; i<=L; ++i) {Nmax[i]=0; idmaxwin[i]=-1;} /* MR1 */
1082 for (k=0; k<N_in; k++) seqid_prev[k]=-1;
1083 if (Ndiff<=0 || Ndiff>=N_in) {seqid1=seqid2; Ndiff=N_in; diffNmax=Ndiff;}
1085 // Check coverage and sim-to-query criteria for each sequence k
1086 for (k=0; k<N_in; k++)
1088 if (keep[k]==KEEP_NOT || keep[k]==KEEP_ALWAYS) continue; // seq k not regular sequence OR is marked sequence
1089 if (100*nres[k]<coverage*L) {keep[k]=KEEP_NOT; continue;} // coverage too low? => reject once and for all
1093 // Check if score-per-column with query is at least qsc
1096 float qsc_min = qsc*nres[k]; // minimum total score of seq k with query
1098 int gapq=0, gapk=0; // number of consecutive gaps in query or k'th sequence at position i
1099 for (int i=first[k]; i<=last[k]; i++)
1104 if (X[kfirst][i]<20)
1107 qsc_sum += S[(int)X[kfirst][i]][(int)X[k][i]];
1109 else if (gapq++) qsc_sum-=PLTY_GAPEXTD; else qsc_sum-=PLTY_GAPOPEN;
1111 else if (X[kfirst][i]<20)
1114 if (gapk++) qsc_sum-=PLTY_GAPEXTD; else qsc_sum-=PLTY_GAPOPEN;
1117 // printf("k=%3i qsc=%6.2f\n",k,qsc_sum);
1118 if (qsc_sum<qsc_min) {keep[k]=KEEP_NOT; continue;} // too different from query? => reject once and for all
1121 //Check if sequence similarity with query at least qid?
1122 if (qdiff_max_frac<0.999)
1124 qdiff_max=int(qdiff_max_frac*nres[k]+0.9999);
1125 // printf("k=%-4i nres=%-4i qdiff_max=%-4i first=%-4i last=%-4i",k,nres[k],qdiff_max,first[k],last[k]);
1127 for (int i=first[k]; i<=last[k]; i++)
1128 // enough different residues to reject based on minimum qid with query? => break
1129 if (X[k][i]<20 && X[k][i]!=X[kfirst][i] && ++diff>=qdiff_max) break;
1130 // printf(" diff=%4i\n",diff);
1131 if (diff>=qdiff_max) {keep[k]=KEEP_NOT; continue;} // too different from query? => reject once and for all
1133 // printf(" qsc=%6.2f qid=%6.2f \n",qsc_sum/nres[k],100.0*(1.0-(float)(diff)/nres[k]));
1138 for (n=k=0; k<N_in; k++) if (keep[k]>KEEP_NOT) n++;
1142 // Successively increment idmax[i] at positons where N[i]<Ndiff
1143 //for (seqid=seqid1; seqid<=seqid2; seqid+=1+(seqid>=50)) /* MR1 */
1145 while (seqid<=seqid2)
1149 for (i=1; i<=L; i++) if (N[i]<Ndiff) idmax[i]=seqid;
1151 // Update idmaxwin[i] as minimum of idmax[i-WFIL,i+WFIL]. If idmaxwin[] has not changed then stop
1153 for (i=1; i<=L; i++)
1155 int idmax_min=seqid2;
1156 for (j=imax(1,imin(L-2*WFIL+1,i-WFIL)); j<=imin(L,imax(2*WFIL,i+WFIL)); j++)
1157 if (idmax[j]<idmax_min) idmax_min=idmax[j];
1158 if (idmax_min>idmaxwin[i]) stop=0; // idmaxwin[i] has changed => do not stop
1159 idmaxwin[i]=idmax_min;
1164 diffNmax_prev = diffNmax;
1166 for (i=1; i<=L; ++i)
1169 for (j=imax(1,imin(L-2*WFIL+1,i-WFIL)); j<=imin(L,imax(2*WFIL,i+WFIL)); ++j)
1170 if (N[j]>max) max=N[j];
1171 if (Nmax[i]<max) Nmax[i]=max;
1176 if (diffNmax<Ndiff-Nmax[i]) diffNmax=Ndiff-Nmax[i];
1181 //printf("seqid=%3i diffNmax_prev= %-4i diffNmax= %-4i n=%-5i N_in-N_ss=%-5i\n",seqid,diffNmax_prev,diffNmax,n,N_in-N_ss);
1186 // printf("idmax ");
1187 // for (i=1; i<=L; i++) printf("%2i ",idmax[i]);
1189 // printf("idmaxwin ");
1190 // for (i=1; i<=L; i++) printf("%2i ",idmaxwin[i]);
1193 // for (i=1; i<=L; i++) printf("%2i ",N[i]);
1196 // Loop over all candidate sequences kk (-> k)
1197 for (kk=0; kk<N_in; kk++)
1199 if (inkk[kk]) continue; // seq k already accepted
1201 if (!keep[k]) continue; // seq k is not regular aa sequence or already suppressed by coverage or qid criterion
1202 if (keep[k]==KEEP_ALWAYS) {inkk[kk]=2; continue;} // accept all marked sequences (no n++, since this has been done already)
1204 // Calculate max-seq-id threshold seqidk for sequence k (as maximum over idmaxwin[i])
1205 if (seqid>=100) {in[k]=inkk[kk]=1; n++; continue;}
1206 float seqidk=seqid1;
1207 for (i=first[k]; i<=last[k]; i++)
1208 if (idmaxwin[i]>seqidk) seqidk=idmaxwin[i];
1209 if (seqid==seqid_prev[k]) continue; // sequence has already been rejected at this seqid threshold => reject this time
1210 seqid_prev[k]=seqid;
1211 diff_min_frac =0.9999-0.01*seqidk; // min fraction of differing positions between sequence j and k needed to accept sequence k
1213 // Loop over already accepted sequences
1214 for (jj=0; jj<kk; jj++)
1216 if (!inkk[jj]) continue;
1218 first_kj=imax(first[k],first[j]);
1219 last_kj =imin(last[k],last[j]);
1220 cov_kj = last_kj-first_kj+1;
1221 diff_suff=int(diff_min_frac*imin(nres[k],cov_kj)+0.999); // nres[j]>nres[k] anyway because of sorting /* MR1 0.999 */
1223 for (int i=first_kj; i<=last_kj; i++)
1225 // enough different residues to accept? => break
1226 if (X[k][i]>=NAA || X[j][i]>=NAA)
1229 if (X[k][i]!=X[j][i] && ++diff>=diff_suff) break; // accept (k,j)
1232 // printf("%20.20s with %20.20s: diff=%i diff_min_frac*cov_kj=%f diff_suff=%i nres=%i cov_kj=%i\n",sname[k],sname[j],diff,diff_min_frac*cov_kj,diff_suff,nres[k],cov_kj);
1233 // printf("%s\n%s\n\n",seq[k],seq[j]);
1235 //if (float(diff)<fmin(diff_min_frac*cov_kj,diff_suff)) break; //similarity > acceptace threshold? Reject! /* MR1 */
1236 if (diff<diff_suff && float(diff)<=diff_min_frac*cov_kj) break; //dissimilarity < acceptace threshold? Reject! /* MR1 */
1240 if (jj>=kk) // did loop reach end? => accept k. Otherwise reject k (the shorter of the two)
1244 for (i=first[k]; i<=last[k]; i++) N[i]++; // update number of sequences at position i
1245 // printf("%i %20.20s accepted\n",k,sname[k]);
1249 // printf("%20.20s rejected: too similar with seq %20.20s diff=%i diff_min_frac*cov_kj=%f diff_suff=%i nres=%i cov_kj=%i\n",sname[k],sname[j],diff,diff_min_frac*cov_kj,diff_suff,nres[k],cov_kj);
1250 // printf("%s\n%s\n\n",seq[k],seq[j]);
1253 } // End Loop over all candidate sequences kk
1257 // printf("seqid_prev[k]= \n");
1258 // for (k=0; k<N_in; k++) printf("%2i ",seqid_prev[k]);
1261 // Increment seqid /* MR1 */
1262 seqid_step = imax(1,imin(5,diffNmax/(diffNmax_prev-diffNmax+1)*seqid_step/2));
1263 seqid += seqid_step;
1265 } // End Loop over seqid
1269 printf("%i out of %i sequences passed filter (",n,N_in-N_ss);
1271 printf("%i%% min coverage, ",coverage);
1273 printf("%i%% min sequence identity to query, ",qid);
1275 printf("%.2f bits min score per column to query, ",qsc);
1276 if (Ndiff<N_in && Ndiff>0)
1277 printf("up to %i%% position-dependent max pairwise sequence identity)\n",seqid);
1279 printf("%i%% max pairwise sequence identity)\n",seqid1);
1282 for (k=0; k<N_in; k++) keep[k]=in[k];
1283 delete[] in; in = NULL;
1284 delete[] inkk; inkk = NULL;
1285 //delete[] idmax; idmax = NULL;
1286 delete[] Nmax; /* MR1 */
1287 delete[] idmaxwin; idmaxwin = NULL;
1288 delete[] seqid_prev; seqid_prev = NULL;
1289 delete[] N; N = NULL;
1291 printf("%s:%s:%d: sequences accepted = %d/%d\n", __FUNCTION__, __FILE__, __LINE__, n, N_in-N_ss);
1298 /* MR1: the Alignment::HomologyFilter is no longer needed in hhalign-stand-alone */
1299 /////////////////////////////////////////////////////////////////////////////
1301 * @brief Filter for min score per column coresc with core query profile,
1302 * defined by coverage_core and qsc_core
1304 /////////////////////////////////////////////////////////////////////////////
1306 Alignment::HomologyFilter(int coverage_core, float qsc_core, float coresc)
1308 const int seqid_core=90; //maximum sequence identity in core alignment
1309 const int qid_core=0;
1310 const int Ndiff_core=0;
1313 char* coreseq=new(char[N_in]); // coreseq[k]=1 if sequence belongs to core of alignment (i.e. it is very similar to query)
1314 for (int k=0; k<N_in; k++) coreseq[k]=keep[k]; // Copy keep[] into coreseq[]
1316 // Remove sequences with seq. identity larger than seqid percent (remove the shorter of two)
1318 n = Filter2(coreseq,coverage_core,qid_core,qsc_core,seqid_core,seqid_core,Ndiff_core);
1322 printf("%i out of %i core alignment sequences passed filter (",n,N_in-N_ss);
1323 if (par.coverage_core)
1324 printf("%i%% min coverage, ",coverage_core);
1326 printf("%i%% min sequence identity to query, ",qid_core);
1328 printf("%.2f bits min score per column to query, ",qsc_core);
1329 printf("%i%% max pairwise sequence identity)\n",seqid_core);
1332 // Calculate bare AA frequencies and transition probabilities -> qcore.f[i][a], qcore.tr[i][a]
1333 FrequenciesAndTransitions(qcore,coreseq);
1335 // Add transition pseudocounts to query -> q.p[i][a] (gapd=1.0, gape=0.333, gapf=gapg=1.0, gaph=gapi=1.0, gapb=1.0
1336 qcore.AddTransitionPseudocounts(1.0,0.333,1.0,1.0,1.0,1.0,1.0);
1338 // Generate an amino acid frequency matrix from f[i][a] with full pseudocount admixture (tau=1) -> g[i][a]
1339 qcore.PreparePseudocounts();
1341 // Add amino acid pseudocounts to query: qcore.p[i][a] = (1-tau)*f[i][a] + tau*g[i][a]
1342 qcore.AddAminoAcidPseudocounts(2,1.5,2.0,1.0); // pcm=2, pca=1.0, pcb=2.5, pcc=1.0
1344 // Filter out all sequences below min score per column with qcore
1345 n=FilterWithCoreHMM(keep, coresc, qcore);
1347 if (v>=2) cout<<n<<" out of "<<N_in-N_ss<<" sequences filtered by minimum score-per-column threshold of "<<qsc_core<<"\n";
1348 delete[] coreseq; coreseq = NULL;
1353 /////////////////////////////////////////////////////////////////////////////////////
1355 * @brief Filter out all sequences below a minimum score per column with profile qcore
1358 Alignment::FilterWithCoreHMM(char in[], float coresc, HMM& qcore)
1360 int k; // count sequences in alignment
1361 int i; // column in query alignment
1362 int a; // amino acid (0..19)
1363 int n=1; // number of sequences that passed filter
1364 float** logodds=new(float*[L+1]); // log-odds ratios for HMM qcore
1365 char gap; // 1: previous state in seq k was a gap 0: previous state in seq k was an amino acid
1366 float score; // score of sequence k aligned with qcore
1368 for (i=1; i<=L; i++) logodds[i]=new(float[21]);
1370 // Determine first[k], last[k]?
1373 first=new(int[N_in]);// first non-gap position in sequence k
1374 last =new(int[N_in]);// last non-gap position in sequence k
1375 for (k=0; k<N_in; k++) // do this for ALL sequences, not only those with in[k]==1 (since in[k] may be display[k])
1377 for (i=1; i<=L; i++) if (X[k][i]<NAA) break;
1379 for (i=L; i>=1; i--) if (X[k][i]<NAA) break;
1384 // Determine number of residues nres[k]?
1387 nres=new(int[N_in]);
1388 for (k=0; k<N_in; k++) // do this for ALL sequences, not only those with in[k]==1 (since in[k] may be display[k])
1391 for (i=first[k]; i<=last[k]; i++)
1392 if (X[k][i]<NAA) nr++;
1394 // printf("%20.20s nres=%3i first=%3i last=%3i\n",sname[k],nr,f,l);
1398 // Precalculate the log-odds for qcore
1399 for (i=1; i<=L; i++)
1401 for (a=0; a<NAA; a++)
1402 logodds[i][a]=fast_log2(qcore.p[i][a]/pb[a]);
1403 logodds[i][ANY]=-0.5; // half a bit penalty for X
1405 // printf(" A R N D C Q E G H I L K M F P S T W Y V\n");
1406 // printf("%6i ",i);
1407 // for (a=0; a<20; ++a) fprintf(stdout,"%5.1f ",100*qcore.f[i][a]);
1410 // for (a=0; a<20; ++a) fprintf(stdout,"%5.1f ",100*qcore.g[i][a]);
1413 // for (a=0; a<20; ++a) fprintf(stdout,"%5.1f ",100*qcore.p[i][a]);
1416 // for (a=0; a<20; ++a) fprintf(stdout,"%5.1f ",100*pb[a]);
1419 // for (a=0; a<20; ++a) fprintf(stdout,"%5.2f ",fast_log2(qcore.p[i][a]/pb[a]));
1423 // Main loop: test all sequences k
1424 for (k=kfirst+1; k<N_in; k++)
1426 if (!in[k]) continue; // if in[k]==0 sequence k will be suppressed directly
1429 float score_prev=0.0;
1431 // Calculate score of sequence k with core HMM
1433 for (i=first[k]; i<=last[k]; i++)
1436 if (X[k][i]<=ANY) // current state is Match
1438 score_M=logodds[i][ (int)X[k][i]];
1439 score+=logodds[i][ (int)X[k][i]];
1440 if (gap) score+=qcore.tr[i][D2M]; else score+=qcore.tr[i][M2M];
1443 else if (X[k][i]==GAP) // current state is Delete (ignore ENDGAPs)
1445 if (gap) score+=qcore.tr[i][D2D]; else score+=qcore.tr[i][M2D];
1448 if (I[k][i]) score+=qcore.tr[i][M2I]+(I[k][i]-1)*qcore.tr[i][I2I]+qcore.tr[i][I2M];
1449 // if (k==2) printf("i=%3i %c:%c score_M=%6.2f score=%6.2f score_sum=%6.2f \n",i,i2aa(X[kfirst][i]),i2aa(X[k][i]),score_M,score-score_prev,score);
1453 printf("k=%3i score=%6.2f\n",k,score);
1454 if (score<nres[k]*coresc) in[k]=0; else n++;// reject sequence k?
1456 for (i=1; i<=L; i++){
1457 delete[] logodds[i]; logodds[i] = NULL;
1459 delete[] logodds; logodds = NULL;
1466 /////////////////////////////////////////////////////////////////////////////////////
1468 * @brief Filter alignment to given diversity/Neff
1471 Alignment::FilterNeff()
1475 const float TOLX=0.001;
1476 const float TOLY=0.02;
1478 for (int k=0; k<N_in; ++k) dummy[k]=keep[k];
1482 float y0=filter_by_qsc(x0,dummy);
1483 float y1=filter_by_qsc(x1,dummy);
1485 while (y0-par.Neff>0 && par.Neff-y1>0)
1487 x = x0 + (par.Neff-y0)*(x1-x0)/(y1-y0); // linear interpolation between (x0,y0) and (x1,y1)
1488 y = filter_by_qsc(x,dummy);
1489 if (v>=2) printf(" %3i x0=%6.3f -> %6.3f x=%6.3f -> %6.3f x1=%6.3f -> %6.3f \n",++i,x0,y0,x,y,x1,y1);
1490 if (y>par.Neff) {x0=x; y0=y;} else {x1=x; y1=y;}
1491 if (fabs(par.Neff-y)<TOLY || x1-x0<TOLX) break;
1495 if (y0>=par.Neff && y1<=par.Neff)
1497 // Write filtered alignment WITH insert states (lower case) to alignment file
1498 if (v>=2) printf("Found Neff=%6.3f at filter threshold qsc=%6.3f\n",y,x);
1502 printf("Diversity of unfiltered alignment %.2f is below target diversity %.2f. No alignment written\n",y0,par.Neff);
1507 float Alignment::filter_by_qsc(float qsc, char* dummy)
1510 for (int k=0; k<N_in; ++k) keep[k]=dummy[k];
1511 Filter2(keep,par.coverage,0,qsc,par.max_seqid+1,par.max_seqid,0);
1512 FrequenciesAndTransitions(q);
1513 // printf("qsc=%4.1f N_filtered=%-3i Neff=%6.3f\n",qsc,n,q.Neff_HMM);
1518 /////////////////////////////////////////////////////////////////////////////////////
1520 * @brief Calculate AA frequencies q.p[i][a] and transition probabilities q.tr[i][a] from alignment
1523 Alignment::FrequenciesAndTransitions(HMM& q, char* in)
1525 int k; // index of sequence
1526 int i; // position in alignment
1527 int a; // amino acid (0..19)
1528 int ni[NAA+3]; // number of times amino acid a occurs at position i
1529 int naa; // number of different amino acids
1532 cout<<"Calculating position-dependent weights on subalignments\n";
1534 if (in==NULL) in=keep; // what's this good for?
1538 for (k=0; k<N_in; k++) wg[k]=0.0; // initialized wg[k]
1539 // Calculate global weights
1540 for (i=1; i<=L; i++) // for all positions i in alignment
1542 for (a=0; a<20; a++) ni[a]=0;
1543 for (k=0; k<N_in; k++) if (in[k]) ni[ (int)X[k][i]]++;
1544 naa=0; for (a=0; a<20; a++) if(ni[a]) naa++;
1545 if (!naa) naa=1; //naa=0 when column consists of only gaps and Xs (=ANY)
1546 for (k=0; k<N_in; k++)
1547 if (in[k] && X[k][i]<20)
1548 wg[k] += 1.0/float(ni[ (int)X[k][i]]*naa*(nres[k]+30));
1549 // ensure that each residue of a short sequence contributes as much as a residue of a long sequence:
1550 // contribution is proportional to one over sequence length nres[k] plus 30.
1552 NormalizeTo1(wg,N_in);
1555 // Do pos-specific sequence weighting and calculate amino acid frequencies and transitions
1556 for (k=0; k<N_in; k++) X[k][0]=ENDGAP; // make sure that sequences ENTER subalignment j for j=1
1557 for (k=0; k<N_in; k++) X[k][L+1]=ENDGAP; // does it have an influence?
1562 #pragma omp parallel sections
1565 Amino_acid_frequencies_and_transitions_from_M_state(q,in); // use subalignments of seqs with residue in i
1567 Transitions_from_I_state(q,in); // use subalignments of seqs with insert in i
1569 Transitions_from_D_state(q,in); // use subalignments of seqs with delete in i. Must be last of these three calls if par.wg==1!
1574 #pragma omp parallel sections
1577 Amino_acid_frequencies_and_transitions_from_M_state(q,in); // use subalignments of seqs with residue in i
1579 Transitions_from_I_state(q,in); // use subalignments of seqs with insert in i
1581 Transitions_from_D_state(q,in); // use subalignments of seqs with delete in i. Must be last of these three calls if par.wg==1!
1584 Amino_acid_frequencies_and_transitions_from_M_state(q,in);
1585 Transitions_from_I_state(q,in);
1586 Transitions_from_D_state(q,in);
1589 else // N_filtered==1
1591 X[kfirst][0]=X[kfirst][L+1]=ANY; // (to avoid anallowed access within loop)
1593 for (i=0; i<=L+1; i++) // for all positions i in alignment
1596 q.Neff_I[i]=q.Neff_D[i]=0.0f;
1597 for (a=0; a<20; a++) q.f[i][a]=0.0;
1598 /* this is the crucial change that makes terminal-X work */
1599 //q.f[i][ (int)(X[kfirst][i]) ] = 1.0; /* MR1 */
1600 if (X[kfirst][i] < ANY) /* MR1 */
1601 q.f[i][(unsigned int) X[kfirst][i] ] = 1.0;
1603 for (a=0; a<20; ++a) q.f[i][a]=pb[a];
1605 q.tr[i][M2I]=-100000.0;
1606 q.tr[i][M2D]=-100000.0;
1607 q.tr[i][I2M]=-100000.0;
1608 q.tr[i][I2I]=-100000.0;
1609 q.tr[i][D2M]=-100000.0;
1610 q.tr[i][D2D]=-100000.0;
1615 q.Neff_M[0]=q.Neff_I[0]=q.Neff_D[0]=99.999; // Neff_av[0] is used for calculation of transition pseudocounts for the start state
1620 printf("\nMatches:\n");
1621 printf("col Neff nseqs\n");
1622 for (i=1; i<=imin(L,100); i++)
1623 printf("%3i %5.2f %3i\n",i,q.Neff_M[i],nseqs[i]);
1625 printf("\nInserts:\n");
1626 printf("col Neff nseqs\n");
1627 for (i=1; i<=imin(L,100); i++)
1628 printf("%3i %5.2f %3i\n",i,q.Neff_I[i],nseqs[i]);
1630 printf("\nDeletes:\n");
1631 printf("col Neff nseqs\n");
1632 for (i=1; i<=imin(L,100); i++)
1633 printf("%3i %5.2f %3i\n",i,q.Neff_D[i],nseqs[i]);
1636 // Copy column information into HMM q
1639 q.N_filtered=N_filtered;
1640 for (i=1; i<=L; i++) q.l[i]=l[i];
1642 // Set names in HMM q
1643 if (strlen(q.name)==0) strcpy(q.name,name);
1644 if (strlen(q.longname)==0) strcpy(q.longname,longname);
1645 if (strlen(q.fam)==0) strcpy(q.fam,fam);
1646 ScopID(q.cl,q.fold,q.sfam,q.fam); // derive superfamily, fold and class code from family name
1647 strcpy(q.file,file); // Store basename of alignment file name in q.file
1649 // Copy sequences to be displayed into HMM
1650 q.nss_dssp=q.nsa_dssp=q.nss_pred=q.nss_conf=q.nfirst=-1;
1652 if (kss_dssp>=0) q.nss_dssp=n++; // copy dssp sequence?
1653 if (ksa_dssp>=0) q.nsa_dssp=n++; // copy dssp sequence?
1654 if (kss_pred>=0) q.nss_pred=n++; // copy psipred sequence?
1655 if (kss_conf>=0) q.nss_conf=n++; // copy confidence value sequence?
1657 // Calculate consensus sequence?
1658 if (par.showcons || par.cons)
1664 // Reserve space for consensus/conservation sequence as Q-T alignment mark-up
1666 q.sname[q.ncons]=new(char[10]);
1667 if (!q.sname[q.ncons]) {MemoryError("array of names for displayed sequences");}
1668 strcpy(q.sname[q.ncons],"Consensus");
1669 q.seq[q.ncons]=new(char[L+2]);
1670 if (!q.seq[q.ncons]) {MemoryError("array of names for displayed sequences");}
1674 // Reserve space for consensus sequence as first sequence in alignment
1675 q.nfirst=n++; kfirst=-1;
1676 q.sname[q.nfirst]=new(char[strlen(name)+11]);
1677 if (!q.sname[q.nfirst]) {MemoryError("array of names for displayed sequences");}
1678 strcpy(q.sname[q.nfirst],name);
1679 strcat(q.sname[q.nfirst],"_consensus");
1680 q.seq[q.nfirst]=new(char[L+2]);
1681 if (!q.seq[q.nfirst]) {MemoryError("array of names for displayed sequences");}
1683 // Calculate consensus amino acids using similarity matrix
1684 for (i=1; i<=L; i++)
1687 for (a=0; a<20; a++)
1688 if (q.f[i][a]-pb[a]>maxw) {maxw = q.f[i][a]-pb[a]; maxa = a;}
1693 for (int b=0; b<20; b++) maxw += q.f[i][b]*Sim[maxa][b]*Sim[maxa][b];
1694 maxw *= q.Neff_M[i]/(q.Neff_HMM+1); // columns with many gaps don't get consensus symbol
1695 if (maxw>0.6) q.seq[q.ncons][i] = uprchr(i2aa(maxa));
1696 else if (maxw>0.4) q.seq[q.ncons][i] = lwrchr(i2aa(maxa));
1697 else q.seq[q.ncons][i] = 'x';
1699 if (par.cons) q.seq[q.nfirst][i] = uprchr(i2aa(maxa));
1703 q.seq[q.ncons][0]='-';
1704 q.seq[q.ncons][L+1]='\0';
1708 q.seq[q.nfirst][0]='-';
1709 q.seq[q.nfirst][L+1]='\0';
1713 // Copy sequences to be displayed from alignment to HMM
1714 for (k=0; k<N_in; k++)
1719 if (0 && (n>=MAXSEQDIS)) {
1720 /* FIXME: the test was if(n>=MAXSEQDIS),
1721 this test was necessary because alignment memory was static,
1722 now it should be dynamic, and should always have the right size,
1723 there are at least number-of-sequences plus a 'bit' more
1724 however, I do not know what that 'bit' is likely to be (in the future).
1725 at the moment it is 1 for the consnseus and 1 for structure,
1726 but this might change (FS)
1728 if (par.mark) cerr<<"WARNING: maximum number "<<MAXSEQDIS<<" of sequences for display of alignment exceeded\n";
1731 if (k==kss_dssp) nn=q.nss_dssp; // copy dssp sequence to nss_dssp
1732 else if (k==ksa_dssp) nn=q.nsa_dssp;
1733 else if (k==kss_pred) nn=q.nss_pred;
1734 else if (k==kss_conf) nn=q.nss_conf;
1735 else if (k==kfirst) nn=q.nfirst=n++;
1737 // strcut(sname[k]," "); // delete rest of name line beginning with two spaces " " // Why this?? Problem for pdb seqs without chain
1738 q.sname[nn]=new(char[strlen(sname[k])+1]);
1739 if (!q.sname[nn]) {MemoryError("array of names for displayed sequences");}
1740 strcpy(q.sname[nn],sname[k]);
1741 q.seq[nn]=new(char[strlen(seq[k])+1]);
1742 if (!q.seq[nn]) {MemoryError("array of names for displayed sequences");}
1743 strcpy(q.seq[nn],seq[k]);
1746 q.n_display=n; // how many sequences to be displayed in alignments?
1748 // Copy secondary structure information into HMM
1750 for (i=1; i<=L; i++) q.ss_dssp[i]=X[kss_dssp][i];
1752 for (i=1; i<=L; i++) q.sa_dssp[i]=X[ksa_dssp][i];
1755 for (i=1; i<=L; i++) q.ss_pred[i]=X[kss_pred][i];
1757 for (i=1; i<=L; i++) q.ss_conf[i]=X[kss_conf][i];
1759 for (i=1; i<=L; i++) q.ss_conf[i]=5;
1765 // Debug: print occurence of amino acids for each position i
1766 if (v>=2) printf("Effective number of sequences exp(entropy) = %-4.1f\n",q.Neff_HMM); //PRINT
1770 for (a=0; a<20; a++) printf("%4.1f ",100*pb[a]);
1771 cout<<"\nAmino acid frequencies without pseudocounts:\n";
1772 cout<<" A R N D C Q E G H I L K M F P S T W Y V\n";
1773 for (i=1; i<=L; i++)
1776 for (a=0; a<20; a++) printf("%4.0f ",100*q.f[i][a]);
1781 printf("\nListing transition probabilities without pseudocounts:\n");
1782 printf(" i M->M M->I M->D I->M I->I D->M D->D Neff_M Neff_I Neff_D\n");
1783 for (i=0; i<=L; i++)
1785 printf("%4i %6.3f %6.3f %6.3f ",i,pow(2.0,q.tr[i][M2M]),pow(2.0,q.tr[i][M2I]),pow(2.0,q.tr[i][M2D]));
1786 printf("%6.3f %6.3f ",pow(2.0,q.tr[i][I2M]),pow(2.0,q.tr[i][I2I]));
1787 printf("%6.3f %6.3f ",pow(2.0,q.tr[i][D2M]),pow(2.0,q.tr[i][D2D]));
1788 printf("%6.3f %6.3f %6.3f\n",q.Neff_M[i],q.Neff_I[i],q.Neff_D[i]);
1792 q.has_pseudocounts=false; /* MR1 */
1798 /////////////////////////////////////////////////////////////////////////////////////
1800 * FIXME: one of the most time consuming routines (according to gprof on r112)
1803 * @brief Calculate freqs q.f[i][a] and transitions q.tr[i][a] (a=MM,MI,MD) with pos-specific subalignments
1804 * Pos-specific weights are calculated like in "GetPositionSpecificWeights()"
1807 Alignment::Amino_acid_frequencies_and_transitions_from_M_state(HMM& q, char* in)
1809 // Calculate position-dependent weights wi[k] for each i.
1810 // For calculation of weights in column i use sub-alignment
1811 // over sequences which have a *residue* in column i (no gap, no end gap)
1812 // and over columns where none of these sequences has an end gap.
1813 // This is done by updating the arrays n[j][a] at each step i-1->i while letting i run from 1 to L.
1814 // n[j][a] = number of occurences of amino acid a at column j of the subalignment,
1815 // => only columns with n[j][ENDGAP]=0 are contained in the subalignment!
1816 // If no sequences enter or leave the subalignment at the step i-1 -> i (i.e. change=0)
1817 // then the old values wi[k], Neff[i-1], and ncol are used for the new position i.
1818 // Index a can be an amino acid (0-19), ANY=20, GAP=21, or ENDGAP=22
1820 int k; // index of sequence
1821 int i,j; // position in alignment
1822 int a; // amino acid (0..19)
1823 int naa; // number of different amino acids
1824 int** n; // n[j][a] = number of seq's with some residue at column i AND a at position j
1825 //float wi[MAXSEQ]; // weight of sequence k in column i, calculated from subalignment i
1826 float *wi=NULL; // weight of sequence k in column i, calculated from subalignment i
1827 //float Neff[MAXRES]; // diversity of subalignment i
1828 float *Neff = new(float[par.maxResLen]); // diversity of subalignment i
1829 int nseqi=0; // number of sequences in subalignment i
1830 int ncol=0; // number of columns j that contribute to Neff[i]
1831 char change; // has the set of sequences in subalignment changed? 0:no 1:yes
1832 float fj[NAA+3]; // to calculate entropy
1835 wi = new(float[N_in+2]);
1839 for (k=0; k<N_in; k++) wi[k]=wg[k];
1843 Neff[0]=0.0; // if the first column has no residues (i.e. change==0), Neff[i]=Neff[i-1]=Neff[0]
1845 for (j=1; j<=L; j++) n[j]=new(int[NAA+3]);
1846 for (j=1; j<=L; j++)
1847 for (a=0; a<NAA+3; a++) n[j][a]=0;
1850 //////////////////////////////////////////////////////////////////////////////////////////////
1851 // Main loop through alignment columns
1852 for (i=1; i<=L; i++) // Calculate wi[k] at position i as well as Neff[i]
1859 // Check all sequences k and update n[j][a] and ri[j] if necessary
1860 for (k=0; k<N_in; k++)
1862 if (!in[k]) continue;
1863 if (X[k][i-1]>=ANY && X[k][i]<ANY)
1864 { // ... if sequence k was NOT included in i-1 and has to be included for column i
1867 for (int j=1; j<=L; j++) n[j][ (int)X[k][j]]++;
1869 else if (X[k][i-1]<ANY && X[k][i]>=ANY)
1870 { // ... if sequence k WAS included in i-1 and has to be thrown out for column i
1873 for (int j=1; j<=L; j++) n[j][ (int)X[k][j]]--;
1878 // If subalignment changed: update weights wi[k] and Neff[i]
1881 // Initialize weights and numbers of residues for subalignment i
1883 for (k=0; k<N_in; k++) wi[k]=1E-8; // for pathological alignments all wi[k] can get 0; /* MR1 */
1885 // sum wi[k] over all columns j and sequences k of subalignment
1886 for (j=1; j<=L; j++)
1888 // do at least a fraction MAXENDGAPFRAC of sequences in subalignment contain an end gap in j?
1889 if (n[j][ENDGAP]>MAXENDGAPFRAC*nseqi) continue;
1890 naa=0; for (a=0; a<20; a++) if(n[j][a]) naa++;
1891 if (naa==0) continue;
1893 for (k=0; k<N_in; k++)
1895 if (in[k] && X[k][i]<ANY && X[k][j]<ANY)
1897 // if (!n[j][ (int)X[k][j]]) {fprintf(stderr,"Error: Mi=%i: n[%i][X[%i]]=0! (X[%i]=%i)\n",i,j,k,k,X[k][j]);}
1898 wi[k]+=1.0/float(n[j][ (int)X[k][j] ]*naa);
1903 // Check whether number of columns in subalignment is sufficient
1905 // Take global weights
1906 for (k=0; k<N_in; k++)
1907 if(in[k] && X[k][i]<ANY) wi[k]=wg[k]; else wi[k]=0.0;
1909 // Calculate Neff[i]
1911 for (j=1; j<=L; j++)
1913 // do at least a fraction MAXENDGAPFRA of sequences in subalignment contain an end gap in j?
1914 if (n[j][ENDGAP]>MAXENDGAPFRAC*nseqi) continue;
1915 for (a=0; a<20; a++) fj[a]=0;
1916 for (k=0; k<N_in; k++)
1917 if (in[k] && X[k][i]<ANY && X[k][j]<ANY)
1918 fj[ (int)X[k][j] ]+=wi[k];
1919 NormalizeTo1(fj,NAA);
1920 for (a=0; a<20; a++)
1921 if (fj[a]>1E-10) Neff[i]-=fj[a]*fast_log2(fj[a]);
1923 if (ncol>0) Neff[i]=pow(2.0,Neff[i]/ncol); else Neff[i]=1.0;
1927 else //no update was necessary; copy values for i-1
1934 // Calculate amino acid frequencies q.f[i][a] from weights wi[k]
1935 for (a=0; a<20; a++) q.f[i][a]=0;
1936 for (k=0; k<N_in; k++) if (in[k]) q.f[i][ (int)X[k][i] ]+=wi[k];
1937 NormalizeTo1(q.f[i],NAA,pb);
1939 // Calculate transition probabilities from M state
1940 q.tr[i][M2M]=q.tr[i][M2D]=q.tr[i][M2I]=0.0;
1941 for (k=0; k<N_in; k++) //for all sequences
1943 if (!in[k]) continue;
1944 //if input alignment is local ignore transitions from and to end gaps
1945 if (X[k][i]<ANY) //current state is M
1947 if (I[k][i]) //next state is I
1948 q.tr[i][M2I]+=wi[k];
1949 else if (X[k][i+1]<=ANY) //next state is M
1950 q.tr[i][M2M]+=wi[k];
1951 else if (X[k][i+1]==GAP) //next state is D
1952 q.tr[i][M2D]+=wi[k];
1955 // Normalize and take log
1956 sum = q.tr[i][M2M]+q.tr[i][M2I]+q.tr[i][M2D]+FLT_MIN;
1957 q.tr[i][M2M]=log2(q.tr[i][M2M]/sum);
1958 q.tr[i][M2I]=log2(q.tr[i][M2I]/sum);
1959 q.tr[i][M2D]=log2(q.tr[i][M2D]/sum);
1961 // for (k=0; k<N_in; k++) if (in[k]) w[k][i]=wi[k];
1963 // DD TODO:fill in all the missing Neff values
1966 // end loop through alignment columns i
1967 //////////////////////////////////////////////////////////////////////////////////////////////
1969 delete[](wi); wi=NULL;
1971 for (j=1; j<=L; j++){
1972 delete[](n[j]); (n[j]) = NULL;
1974 delete[](n); (n) = NULL;
1977 q.tr[0][M2I]=-100000;
1978 q.tr[0][M2D]=-100000;
1980 q.tr[L][M2I]=-100000;
1981 q.tr[L][M2D]=-100000;
1982 q.Neff_M[0]=99.999; // Neff_av[0] is used for calculation of transition pseudocounts for the start state
1984 // Set emission probabilities of zero'th (begin) state and L+1st (end) state to background probabilities
1985 for (a=0; a<20; a++) q.f[0][a]=q.f[L+1][a]=pb[a];
1987 // Assign Neff_M[i] and calculate average over alignment, Neff_M[0]
1990 for (i=1; i<=L; i++)
1993 for (a=0; a<20; a++)
1994 if (q.f[i][a]>1E-10) sum -= q.f[i][a]*fast_log2(q.f[i][a]);
1995 q.Neff_HMM+=pow(2.0,sum);
1998 float Nlim=fmax(10.0,q.Neff_HMM+1.0); // limiting Neff
1999 float scale=log2((Nlim-q.Neff_HMM)/(Nlim-1.0)); // for calculating Neff for those seqs with inserts at specific pos
2000 for (i=1; i<=L; i++)
2002 float w_M=-1.0/N_filtered;
2003 for (k=0; k<N_in; k++)
2004 if (in[k] && X[k][i]<=ANY) w_M+=wg[k];
2005 if (w_M<0) q.Neff_M[i]=1.0;
2006 else q.Neff_M[i] = Nlim - (Nlim-1.0)*fpow2(scale*w_M);
2007 // fprintf(stderr,"M i=%3i ncol=--- Neff_M=%5.2f Nlim=%5.2f w_M=%5.3f Neff_M=%5.2f\n",i,q.Neff_HMM,Nlim,w_M,q.Neff_M[i]);
2012 for (i=1; i<=L; i++)
2014 q.Neff_HMM+=Neff[i];
2015 q.Neff_M[i]=Neff[i];
2016 if (q.Neff_M[i] == 0) { q.Neff_M[i] = 1; } /* MR1 */
2021 delete[] Neff; Neff = NULL;
2025 } /* this is the end of Alignment::Amino_acid_frequencies_and_transitions_from_M_state() */
2028 /////////////////////////////////////////////////////////////////////////////////////
2030 * @brief Calculate transitions q.tr[i][a] (a=DM,DD) with pos-specific subalignments
2033 Alignment::Transitions_from_I_state(HMM& q, char* in)
2035 // Calculate position-dependent weights wi[k] for each i.
2036 // For calculation of weights in column i use sub-alignment
2037 // over sequences which have a INSERT in column i
2038 // and over columns where none of these sequences has an end gap.
2039 // This is done by calculating the arrays n[j][a] and rj[j] at each step i-1->i while letting i run from 1 to L.
2040 // n[j][a] = number of occurences of amino acid a at column j of the subalignment,
2041 // => only columns with n[j][ENDGAP]=0 are contained in the subalignment!
2042 // If no sequences enter or leave the subalignment at the step i-1 -> i (i.e. change=0)
2043 // then the old values wi[k], Neff[i-1], and ncol are used for the new position i.
2044 // Index a can be an amino acid (0-19), ANY=20, GAP=21, or ENDGAP=22
2046 int k; // index of sequence
2047 int i,j; // position in alignment
2048 int a; // amino acid (0..19)
2049 int naa; // number of different amino acids
2050 int** n; // n[j][a] = number of seq's with some residue at column i AND a at position j
2051 //float wi[MAXSEQ]; // weight of sequence k in column i, calculated from subalignment i
2052 float *wi = NULL; // weight of sequence k in column i, calculated from subalignment i
2053 //float Neff[MAXRES]; // diversity of subalignment i
2054 float *Neff = new(float[par.maxResLen]); // diversity of subalignment i
2055 int nseqi; // number of sequences in subalignment i
2056 int ncol; // number of columns j that contribute to Neff[i]
2057 float fj[NAA+3]; // to calculate entropy
2059 float Nlim=0.0; // only for global weights
2060 float scale=0.0; // only for global weights
2062 wi = new(float[N_in+2]);
2067 for (k=0; k<N_in; k++) wi[k]=wg[k];
2068 Nlim=fmax(10.0,q.Neff_HMM+1.0); // limiting Neff
2069 scale=log2((Nlim-q.Neff_HMM)/(Nlim-1.0)); // for calculating Neff for those seqs with inserts at specific pos
2074 for (j=1; j<=L; j++) n[j]=new(int[NAA+3]);
2076 //////////////////////////////////////////////////////////////////////////////////////////////
2077 // Main loop through alignment columns
2078 for (i=1; i<=L; i++) // Calculate wi[k] at position i as well as Neff[i]
2080 if (par.wg==0) // local weights?
2083 // Calculate n[j][a] and ri[j]
2085 for (k=0; k<N_in; k++)
2087 if (in[k] && I[k][i]>0)
2089 if (nseqi==0) // Initialize only if inserts present! Otherwise O(L*L) even for single sequences!
2091 // Initialization of n[j][a]
2092 for (j=1; j<=L; j++)
2093 for (a=0; a<NAA+3; a++) n[j][a]=0;
2096 for (int j=1; j<=L; j++) n[j][ (int)X[k][j]]++;
2101 // If there is no sequence in subalignment j ...
2105 Neff[i]=0.0; // effective number of sequence = 0!
2106 q.tr[i][I2M]=-100000;
2107 q.tr[i][I2I]=-100000;
2111 // update weights wi[k] and Neff[i]
2114 // Initialize weights and numbers of residues for subalignment i
2116 for (k=0; k<N_in; k++) wi[k]=0.0;
2118 // sum wi[k] over all columns j and sequences k of subalignment
2119 for (j=1; j<=L; j++)
2121 if (n[j][ENDGAP]>MAXENDGAPFRAC*nseqi) continue;
2122 naa=0; for (a=0; a<20; a++) if(n[j][a]) naa++;
2123 if (naa==0) continue;
2125 for (k=0; k<N_in; k++)
2127 if (in[k] && I[k][i]>0 && X[k][j]<ANY)
2129 if (!n[j][ (int)X[k][j]]) {fprintf(stderr,"Error: Ii=%i: n[%i][X[%i]]=0! (X[%i]=%i)\n",i,j,k,k,X[k][j]);}
2130 wi[k]+=1.0/float(n[j][ (int)X[k][j] ]*naa);
2135 // Check whether number of columns in subalignment is sufficient
2137 // Take global weights
2138 for (k=0; k<N_in; k++)
2139 if(in[k] && I[k][i]>0) wi[k]=wg[k]; else wi[k]=0.0;
2141 // Calculate Neff[i]
2143 for (j=1; j<=L; j++)
2145 if (n[j][ENDGAP]>MAXENDGAPFRAC*nseqi) continue;
2146 for (a=0; a<20; a++) fj[a]=0;
2147 for (k=0; k<N_in; k++)
2148 if (in[k] && I[k][i]>0 && X[k][j]<ANY)
2149 fj[ (int)X[k][j] ]+=wi[k];
2150 NormalizeTo1(fj,NAA);
2151 for (a=0; a<20; a++)
2152 if (fj[a]>1E-10) Neff[i]-=fj[a]*fast_log2(fj[a]);
2154 if (ncol>0) Neff[i]=pow(2.0,Neff[i]/ncol); else Neff[i]=1.0;
2157 // Calculate transition probabilities from I state
2158 q.tr[i][I2M]=q.tr[i][I2I]=0.0;
2159 for (k=0; k<N_in; k++) //for all sequences
2161 if (in[k] && I[k][i]>0) //current state is I
2163 q.tr[i][I2M]+=wi[k];
2164 q.tr[i][I2I]+=wi[k]*(I[k][i]-1);
2169 else // fast global weights?
2171 float w_I=-1.0/N_filtered;
2173 q.tr[i][I2M]=q.tr[i][I2I]=0.0;
2174 // Calculate amino acid frequencies fj[a] from weights wg[k]
2175 for (k=0; k<N_in; k++)
2176 if (in[k] && I[k][i]>0)
2180 q.tr[i][I2M]+=wi[k];
2181 q.tr[i][I2I]+=wi[k]*(I[k][i]-1);
2185 if (w_I<0) Neff[i]=1.0;
2186 else Neff[i] = Nlim - (Nlim-1.0)*fpow2(scale*w_I);
2187 // fprintf(stderr,"I i=%3i ncol=%3i Neff_M=%5.2f Nlim=%5.2f w_I=%5.3f Neff_I=%5.2f\n",i,ncol,q.Neff_HMM,Nlim,w_I,Neff[i]);
2192 q.tr[i][I2M]=-100000;
2193 q.tr[i][I2I]=-100000;
2198 // Normalize and take log
2199 sum = q.tr[i][I2M]+q.tr[i][I2I];
2200 q.tr[i][I2M]=log2(q.tr[i][I2M]/sum);
2201 q.tr[i][I2I]=log2(q.tr[i][I2I]/sum);
2204 // end loop through alignment columns i
2205 //////////////////////////////////////////////////////////////////////////////////////////////
2207 delete[](wi); wi = NULL;
2209 for (j=1; j<=L; j++){
2210 delete[](n[j]); (n[j]) = NULL;
2212 delete[](n); (n) = NULL;
2215 q.tr[0][I2I]=-100000;
2217 q.tr[L][I2I]=-100000;
2221 for (i=1; i<=L; i++) // Calculate wi[k] at position i as well as Neff[i] and Neff[i]
2222 q.Neff_I[i]=Neff[i];
2224 delete[] Neff; Neff = NULL;
2227 } /* this is the end of Alignment::Transitions_from_I_state() */
2230 /////////////////////////////////////////////////////////////////////////////////////
2232 * @brief Calculate transitions q.tr[i][a] (a=DM,DD) with pos-specific subalignments
2235 Alignment::Transitions_from_D_state(HMM& q, char* in)
2237 // Calculate position-dependent weights wi[k] for each i.
2238 // For calculation of weights in column i use sub-alignment
2239 // over sequences which have a DELETE in column i
2240 // and over columns where none of these sequences has an end gap.
2241 // This is done by updating the arrays n[j][a] and rj[j] at each step i-1->i while letting i run from 1 to L.
2242 // n[j][a] = number of occurences of index a at column j of the subalignment,
2243 // => only columns with n[j][ENDGAP]=0 are contained in the subalignment!
2244 // If no sequences enter or leave the subalignment at the step i-1 -> i (i.e. change=0)
2245 // then the old values wi[k], Neff[i-1], and ncol are used for the new position i.
2246 // Index a can be an amino acid (0-19), ANY=20, GAP=21, or ENDGAP=22
2248 int k; // index of sequence
2249 int i,j; // position in alignment
2250 int a; // amino acid (0..19)
2251 int naa; // number of different amino acids
2252 int** n; // n[j][a] = number of seq's with some residue at column i AND a at position j
2253 //float wi[MAXSEQ]; // weight of sequence k in column i, calculated from subalignment i
2254 float *wi=NULL; // weight of sequence k in column i, calculated from subalignment i
2255 //float Neff[MAXRES]; // diversity of subalignment i
2256 float *Neff = new(float[par.maxResLen]); // diversity of subalignment i
2257 int nseqi=0; // number of sequences in subalignment i (for DEBUGGING)
2258 int ncol=0; // number of columns j that contribute to Neff[i]
2259 char change; // has the set of sequences in subalignment changed? 0:no 1:yes
2260 float fj[NAA+3]; // to calculate entropy
2262 float Nlim=0.0; // only for global weights
2263 float scale=0.0; // only for global weights
2265 wi = new(float[N_in+2]); /* FIXME: FS */
2269 for (k=0; k<N_in; k++) wi[k]=wg[k];
2270 Nlim=fmax(10.0,q.Neff_HMM+1.0); // limiting Neff
2271 scale=log2((Nlim-q.Neff_HMM)/(Nlim-1.0)); // for calculating Neff for those seqs with dels at specific pos
2276 for (j=1; j<=L; j++) n[j]=new(int[NAA+3]);
2277 for (j=1; j<=L; j++)
2278 for (a=0; a<NAA+3; a++) n[j][a]=0;
2282 //////////////////////////////////////////////////////////////////////////////////////////////
2283 // Main loop through alignment columns
2284 for (i=1; i<=L; i++) // Calculate wi[k] at position i as well as Neff[i]
2286 if (par.wg==0) // if local weights
2289 // Check all sequences k and update n[j][a] and ri[j] if necessary
2290 for (k=0; k<N_in; k++)
2292 if (!in[k]) continue;
2293 if (X[k][i-1]!=GAP && X[k][i]==GAP)
2294 { // ... if sequence k was NOT included in i-1 and has to be included for column i
2297 for (int j=1; j<=L; j++) n[j][ (int)X[k][j]]++;
2299 else if (X[k][i-1]==GAP && X[k][i]!=GAP)
2300 { // ... if sequence k WAS included in i-1 and has to be thrown out for column i
2303 for (int j=1; j<=L; j++) n[j][ (int)X[k][j]]--;
2308 // If there is no sequence in subalignment j ...
2312 Neff[i]=0.0; // effective number of sequences = 0!
2313 q.tr[i][D2M]=-100000;
2314 q.tr[i][D2D]=-100000;
2318 // If subalignment changed: update weights wi[k] and Neff[i]
2321 // Initialize weights and numbers of residues for subalignment i
2323 for (k=0; k<N_in; k++) wi[k]=0.0;
2325 // sum wg[k][i] over all columns j and sequences k of subalignment
2326 for (j=1; j<=L; j++)
2328 if (n[j][ENDGAP]>MAXENDGAPFRAC*nseqi) continue;
2329 naa=0; for (a=0; a<20; a++) if(n[j][a]) naa++;
2330 if (naa==0) continue;
2332 for (k=0; k<N_in; k++)
2334 if (in[k] && X[k][i]==GAP && X[k][j]<ANY)
2336 if (!n[j][ (int)X[k][j]]) {fprintf(stderr,"Error: Di=%i: n[%i][X[%i]]=0! (X[%i]=%i)\n",i,j,k,k,X[k][j]);}
2337 wi[k]+=1.0/float(n[j][ (int)X[k][j] ]*naa);
2342 // Check whether number of columns in subalignment is sufficient
2344 // Take global weights
2345 for (k=0; k<N_in; k++)
2346 if(in[k] && X[k][i]==GAP) wi[k]=wg[k]; else wi[k]=0.0;
2348 // Calculate Neff[i]
2350 for (j=1; j<=L; j++)
2352 if (n[j][ENDGAP]>MAXENDGAPFRAC*nseqi) continue;
2353 for (a=0; a<20; a++) fj[a]=0;
2354 for (k=0; k<N_in; k++)
2355 if (in[k] && X[k][i]==GAP && X[k][j]<ANY)
2356 fj[ (int)X[k][j] ]+=wi[k];
2357 NormalizeTo1(fj,NAA);
2358 for (a=0; a<20; a++)
2359 if (fj[a]>1E-10) Neff[i]-=fj[a]*fast_log2(fj[a]);
2361 if (ncol>0) Neff[i]=pow(2.0,Neff[i]/ncol); else Neff[i]=1.0;
2365 else //no update was necessary; copy values for i-1
2370 // Calculate transition probabilities from D state
2371 q.tr[i][D2M]=q.tr[i][D2D]=0.0;
2372 for (k=0; k<N_in; k++) //for all sequences
2374 if (in[k] && X[k][i]==GAP) //current state is D
2376 if (X[k][i+1]==GAP) //next state is D
2377 q.tr[i][D2D]+=wi[k];
2378 else if (X[k][i+1]<=ANY) //next state is M
2379 q.tr[i][D2M]+=wi[k];
2384 else // fast global weights?
2386 float w_D=-1.0/N_filtered;
2388 q.tr[i][D2M]=q.tr[i][D2D]=0.0;
2389 // Calculate amino acid frequencies fj[a] from weights wg[k]
2390 for (k=0; k<N_in; k++) //for all sequences
2391 if (in[k] && X[k][i]==GAP) //current state is D
2395 if (X[k][i+1]==GAP) //next state is D
2396 q.tr[i][D2D]+=wi[k];
2397 else if (X[k][i+1]<=ANY) //next state is M
2398 q.tr[i][D2M]+=wi[k];
2402 if (w_D<0) Neff[i]=1.0;
2403 else Neff[i] = Nlim - (Nlim-1.0)*fpow2(scale*w_D);
2404 // fprintf(stderr,"D i=%3i ncol=%3i Neff_M=%5.2f Nlim=%5.2f w_D=%5.3f Neff_D=%5.2f\n",i,ncol,q.Neff_HMM,Nlim,w_D,Neff[i]);
2408 Neff[i]=0.0; // effective number of sequences = 0!
2409 q.tr[i][D2M]=-100000;
2410 q.tr[i][D2D]=-100000;
2415 // Normalize and take log
2416 sum = q.tr[i][D2M]+q.tr[i][D2D];
2417 q.tr[i][D2M]=log2(q.tr[i][D2M]/sum);
2418 q.tr[i][D2D]=log2(q.tr[i][D2D]/sum);
2421 // end loop through alignment columns i
2422 //////////////////////////////////////////////////////////////////////////////////////////////
2425 q.tr[0][D2D]=-100000;
2429 for (i=1; i<=L; i++)
2430 q.Neff_D[i]=Neff[i];
2432 delete[](wi); wi = NULL;/* FIXME: FS */
2434 for (j=1; j<=L; j++){
2435 delete[](n[j]); (n[j]) = NULL;
2437 delete[](n); (n) = NULL;
2439 delete[] Neff; Neff = NULL;
2442 } /* this is the end of Alignment::Transitions_from_D_state() */
2446 /////////////////////////////////////////////////////////////////////////////////////
2448 * @brief Write alignment without insert states (lower case) to alignment file?
2451 Alignment::WriteWithoutInsertsToFile(char* alnfile)
2453 if (v>=2) cout<<"Writing alignment to "<<alnfile<<"\n";
2455 if (!par.append) alnf = fopen(alnfile,"w"); else alnf = fopen(alnfile,"a");
2456 if (!alnf) OpenFileError(alnfile);
2457 // If alignment name is different from that of query: write name into commentary line
2458 if (strncmp(longname,sname[kfirst],DESCLEN-1)) fprintf(alnf,"#%s\n",longname);
2459 if (v>=2) cout<<"Writing alignment to "<<alnfile<<"\n";
2460 for (int k=0; k<N_in; k++)
2461 if (keep[k] || display[k]==KEEP_ALWAYS) // print if either in profile (keep[k]>0) or display is obligatory (display[k]==2)
2463 fprintf(alnf,">%s\n",sname[k]);
2464 for (int i=1; i<=L; i++) fprintf(alnf,"%c",i2aa(X[k][i]));
2470 /////////////////////////////////////////////////////////////////////////////////////
2471 // Write stored,filtered sequences WITH insert states (lower case) to alignment file?
2472 /////////////////////////////////////////////////////////////////////////////////////
2473 void Alignment::WriteToFile(char* alnfile, const char format[])
2476 if (!par.append) alnf = fopen(alnfile,"w"); else alnf = fopen(alnfile,"a");
2477 if (!alnf) OpenFileError(alnfile);
2478 // If alignment name is different from that of query: write name into commentary line
2479 if (strncmp(longname,sname[kfirst],DESCLEN-1)) fprintf(alnf,"#%s\n",longname);
2480 if (!format || !strcmp(format,"a3m"))
2482 if (v>=2) cout<<"Writing A3M alignment to "<<alnfile<<"\n";
2483 for (int k=0; k<N_in; k++)
2484 if (keep[k] || display[k]==KEEP_ALWAYS) // print if either in profile (keep[k]>0) or display obligatory (display[k]==2)
2485 fprintf(alnf,">%s\n%s\n",sname[k],seq[k]+1);
2487 else // PSI-BLAST format
2489 if (v>=2) cout<<"Writing PSI-BLAST-formatted alignment to "<<alnfile<<"\n";
2490 for (int k=kfirst; k<N_in; k++) // skip sequences before kfirst!!
2491 if (keep[k] || display[k]==KEEP_ALWAYS) // print if either in profile (keep[k]>0) or display obligatory (display[k]==2)
2494 fprintf(alnf,"%-20.20s ",sname[k]);
2495 // for (int i=1; i<=L; i++) fprintf(alnf,"%c",i2aa(X[k][i]));
2496 // fprintf(alnf,"\n");
2498 for (; *ptr!='\0'; ptr++)
2499 if (*ptr==45 || (*ptr>=65 && *ptr<=90)) fprintf(alnf,"%c",*ptr);
2510 * FIXME: this function contains a reference to MAXSEQ & MAXCOL
2511 * however, this may not be accessed (FS)
2513 /////////////////////////////////////////////////////////////////////////////////////
2515 * @brief Read a3m slave alignment of hit from file and merge into (query) master alignment
2518 Alignment::MergeMasterSlave(Hit& hit, char ta3mfile[])
2521 char* cur_seq = new(char[MAXCOL]); // Sequence currently read in
2523 int l,ll; // position in unaligned template (T) sequence Tali.seq[l]
2524 int i; // counts match states in query (Q) HMM
2525 int j; // counts match states in T sequence Tali.seq[l]
2526 int h; // position in aligned T sequence cur_seq[h]
2527 int k; // sequence index
2529 printf("****************%s:%s:%d: did get into MergeMasterSlave\n", __FUNCTION__, __FILE__, __LINE__);
2530 if (v>=3) printf("Merging %s to query alignment\n",ta3mfile);
2532 // If par.append==1 do not print query alignment
2533 if (par.append) for (k=0; k<N_in; k++) keep[k]=display[k]=KEEP_NOT;
2535 // Read template alignment into Tali
2536 FILE* ta3mf=fopen(ta3mfile,"r");
2537 if (!ta3mf) OpenFileError(ta3mfile);
2538 Tali.Read(ta3mf,ta3mfile);
2541 // Filter Tali alignment
2542 Tali.Compress(ta3mfile);
2543 N_filtered = Tali.Filter(par.max_seqid,par.coverage,par.qid,par.qsc,par.Ndiff);
2546 int* imatch=new(int[hit.j2+1]);
2547 int step = hit.nsteps;
2548 for (j=hit.j1; j<=hit.j2; j++)
2550 // Advance to position of next T match state j
2551 while (hit.j[step]<j) step--;
2552 imatch[j] = hit.i[step];
2553 // printf("step=%-3i i=%-3i j=%-3i\n",step,imatch[j],j);
2556 // Determine number of match states of Qali
2557 for (L=0,l=1; seq[kfirst][l]>'\0'; l++)
2558 if ((seq[kfirst][l]>='A' && seq[kfirst][l]<='Z') || seq[kfirst][l]=='-') L++;
2560 // For each sequence in T alignment: align to Qali
2561 for (k=0; k<Tali.N_in; k++)
2563 if (!Tali.keep[k]) continue;
2566 fprintf(stderr,"WARNING in %s: maximum number of %i sequences exceeded while reading %s. Skipping all following sequences\n",program_name,MAXSEQ,ta3mfile);
2569 cur_seq[0]=' '; // 0'th position not used
2571 // Add the hit.i1-1 left end gaps to aligned sequence
2572 for (h=1; h<hit.i1; h++) cur_seq[h]='-';
2574 // Advance to match state hit.j1 of Tali.seq[k]
2575 for (j=0, l=1; (c=Tali.seq[k][l])>'\0'; l++)
2576 if ((c>='A' && c<='Z') || c=='-') // match state at position l?
2577 if ((++j)==hit.j1) break; // yes: increment j. Reached hit,j1? yes: break
2580 {printf("Error: did not find %i match states in sequence %i of %s. Sequence:\n%s\n",hit.j1,k,Tali.name,Tali.seq[k]); exit(1);}
2582 // Write first match state to cur_seq
2583 int iprev=hit.i1; // index of previous query match state
2584 int lprev=l; // previous T match state in Tali.seq[k][l]
2585 cur_seq[h++] = Tali.seq[k][l]; // first column of alignment is Match-Match state
2587 // For each further match state j in alignment
2589 for (j=hit.j1+1; j<=hit.j2; j++)
2591 // Advance to position of next T match state j
2594 // Advance to position of next T match state j
2595 while ((c=Tali.seq[k][++l])>'\0' && ((c>='a' && c<='z') || c=='.')) ;
2597 int di=i-iprev; // number of Match states in Q between T match state j-1 and j
2598 int dl=l-lprev; // 1 + number of inserted residues in T sequence between T match state j-1 and j
2601 // One Q match state for one T match state (treated as special case for speed reasons)
2603 // Q: XXXXXX.....XXXXXX
2604 // T: YYYYYYyyyyyYYYYYY
2608 // Inserts in lower case
2609 for (ll=lprev+1; ll<l; ll++)
2610 if (Tali.seq[k][ll]!='-' && Tali.seq[k][ll]!='.') cur_seq[h++] = lwrchr(Tali.seq[k][ll]);
2612 // Template Match state -> upper case
2613 cur_seq[h++] = Tali.seq[k][ll];
2617 // Gap in query: no Q match state for on T match state (special case for speed reasons)
2619 // Q: XXXXXX.....~~~XXX
2620 // T: YYYYYYyyyyyYYYYYY
2624 // All T residues (including T match state) in lower case
2625 for (ll=lprev+1; ll<=l; ll++)
2626 if (Tali.seq[k][ll]!='-' && Tali.seq[k][ll]!='.') cur_seq[h++] = lwrchr(Tali.seq[k][ll]);
2630 // More Match states in Q than Inserts in the T sequence
2631 // => half T inserts y left, half right-aligned in uc, gaps to fill up
2632 // Number of T insert residues to be left-aligned: (int)(dl/2)
2634 // Q: XXXXXXXXXXXXXXXXXX
2635 // T: YYYYYYYyyy-yyYYYYY
2639 // Add left-bounded template residues
2640 for (ll=lprev+1; ll<=lprev+(int)(dl/2); ll++)
2641 cur_seq[h++]=uprchr(Tali.seq[k][ll]);
2644 for (int gap=1; gap<=di-dl; gap++) cur_seq[h++]='-';
2646 // Add right-bounded residues
2648 cur_seq[h++]=uprchr(Tali.seq[k][ll]);
2652 // Fewer Match states in Q than inserts in T sequence
2653 // => half of available space di for left- half for right-aligned T inserts, rest in lc
2654 // number of T inserts to be left-aligned in uc: (int)(di/2),
2656 // Q: XXXXXXXXX.XXXXXXX
2657 // T: YYYYYYYyyyyyYYYYY
2661 // Add left-bounded template residues
2662 for (ll=lprev+1; ll<=lprev+(int)(di/2); ll++)
2663 cur_seq[h++]=uprchr(Tali.seq[k][ll]);
2665 // Add central inserts
2666 for (int ins=1; ins<=dl-di; ins++,ll++)
2667 if (Tali.seq[k][ll]!='-' && Tali.seq[k][ll]!='.') cur_seq[h++] = lwrchr(Tali.seq[k][ll]);
2669 // Add right-bounded residues
2671 cur_seq[h++]=uprchr(Tali.seq[k][ll]);
2673 // printf("i=%-3i j=%-3i l=%-3i cur_seq=%s\n",i,j,l,cur_seq);
2676 if (h>=maxcol-1000) // too few columns? Reserve double space
2678 char* new_seq=new(char[2*maxcol]);
2679 strncpy(new_seq,cur_seq,maxcol); //////// check: maxcol-1 ????
2680 delete[](cur_seq); (cur_seq) = NULL;
2686 // Add the remaining gaps '-' to the end of the template sequence
2687 for (i=hit.i2+1; i<=L; i++) cur_seq[h++]='-';
2690 keep[N_in] = display[N_in] = KEEP_CONDITIONALLY;
2691 seq[N_in]=new(char[h]);
2692 if (!seq[N_in]) MemoryError("array for input sequences");
2693 strcpy(seq[N_in],cur_seq);
2694 X[N_in]=new(char[h]);
2695 if (!X[N_in]) MemoryError("array for input sequences");
2696 I[N_in]=new(short unsigned int[h]);
2697 if (!I[N_in]) MemoryError("array for input sequences");
2698 sname[N_in]=new(char[strlen(Tali.sname[k])+1]);
2699 if (!sname[N_in]) MemoryError("array for input sequences");
2700 strcpy(sname[N_in],Tali.sname[k]);
2703 // printf("k=%-3i %s\n",k,Tali.seq[k]);
2704 // printf("Query %s\n",seq[kfirst]);
2705 // printf("k=%-3i %s\n\n",k,cur_seq);
2709 // printf("N_in=%-5i HMM=%s with %i sequences\n",N_in,ta3mfile,N_filtered);
2711 delete[] cur_seq; cur_seq = NULL;
2712 delete[] imatch; imatch = NULL;
2713 delete[] ksort; ksort=NULL; // if ksort already existed it will be to short for merged alignment
2714 delete[] first; first=NULL; // if first already existed it will be to short for merged alignment
2715 delete[] last; last=NULL; // if last already existed it will be to short for merged alignment
2717 } /* this is the end of Alignment::MergeMasterSlave() */
2720 /////////////////////////////////////////////////////////////////////////////////////
2722 * @brief Add a sequence to Qali
2725 Alignment::AddSequence(char Xk[], int Ik[])
2727 int i; // position in query and target
2728 if (L<=0) InternalError("L is not set in AddSequence()");
2729 X[N_in]=new(char[L+2]);
2730 for (i=0; i<=L+1; i++) X[N_in][i]=Xk[i];
2732 for (i=0; i<=L+1; i++) I[N_in][i]=0;
2734 for (i=0; i<=L+1; i++) I[N_in][i]=Ik[i];
2739 /////////////////////////////////////////////////////////////////////////////////////
2741 * @brief Determine matrix of position-specific weights w[k][i] for multiple alignment
2742 * Pos-specific weights are calculated like in "Amino_acid_frequencies_and_transitions_from_M_state()"
2745 Alignment::GetPositionSpecificWeights(float* w[])
2747 // Calculate position-dependent weights wi[k] for each i.
2748 // For calculation of weights in column i use sub-alignment
2749 // over sequences which have a *residue* in column i (no gap, no end gap)
2750 // and over columns where none of these sequences has an end gap.
2751 // This is done by updating the arrays n[j][a] at each step i-1->i while letting i run from 1 to L.
2752 // n[j][a] = number of occurences of amino acid a at column j of the subalignment,
2753 // => only columns with n[j][ENDGAP]=0 are contained in the subalignment!
2754 // If no sequences enter or leave the subalignment at the step i-1 -> i (i.e. change=0)
2755 // then the old values w[k][i] and ncol are used for the new position i.
2756 // Index a can be an amino acid (0-19), ANY=20, GAP=21, or ENDGAP=22
2758 char* in=keep; // to keep the code similar to Amino_acid_frequencies_and_transitions_from_M_state()
2759 int k; // index of sequence
2760 int i,j; // position in alignment
2761 int a; // amino acid (0..19)
2762 int naa; // number of different amino acids
2763 int** n; // n[j][a] = number of seq's with some residue at column i AND a at position j
2764 int nseqi=0; // number of sequences in subalignment i
2765 int ncol=0; // number of columns j that contribute to Neff[i]
2766 char change; // has the set of sequences in subalignment changed? 0:no 1:yes
2772 for (k=0; k<N_in; k++)
2773 for (i=1; i<=L; i++) w[k][i]=wg[k];
2780 for (j=1; j<=L; j++) n[j]=new(int[NAA+3]);
2781 for (j=1; j<=L; j++)
2782 for (a=0; a<NAA+3; a++) n[j][a]=0;
2784 //////////////////////////////////////////////////////////////////////////////////////////////
2785 // Main loop through alignment columns
2786 for (i=1; i<=L; i++) // Calculate w[k][i]
2789 // Check all sequences k and update n[j][a] and ri[j] if necessary
2790 for (k=0; k<N_in; k++)
2792 if (!in[k]) continue;
2793 if (X[k][i-1]>=ANY && X[k][i]<ANY)
2794 { // ... if sequence k was NOT included in i-1 and has to be included for column i
2797 for (int j=1; j<=L; j++) n[j][ (int)X[k][j]]++;
2799 else if (X[k][i-1]<ANY && X[k][i]>=ANY)
2800 { // ... if sequence k WAS included in i-1 and has to be thrown out for column i
2803 for (int j=1; j<=L; j++) n[j][ (int)X[k][j]]--;
2808 // If subalignment changed: update weights w[k][i] and Neff[i]
2811 // Initialize weights and numbers of residues for subalignment i
2813 for (k=0; k<N_in; k++) w[k][i]=0.0;
2815 // sum wi[k] over all columns j and sequences k of subalignment
2816 for (j=1; j<=L; j++)
2818 // do at least a fraction MAXENDGAPFRAC of sequences in subalignment contain an end gap in j?
2819 if (n[j][ENDGAP]>MAXENDGAPFRAC*nseqi) continue;
2820 naa=0; for (a=0; a<20; a++) if(n[j][a]) naa++;
2821 if (naa==0) continue;
2823 for (k=0; k<N_in; k++)
2825 if (in[k] && X[k][i]<ANY && X[k][j]<ANY)
2827 // if (!n[j][ (int)X[k][j]]) {fprintf(stderr,"Error: Mi=%i: n[%i][X[%i]]=0! (X[%i]=%i)\n",i,j,k,k,X[k][j]);}
2828 w[k][i]+=1.0/float(n[j][ (int)X[k][j] ]*naa);
2833 // Check whether number of columns in subalignment is sufficient
2835 // Take global weights
2836 for (k=0; k<N_in; k++)
2837 if(in[k]) {if(X[k][i]<ANY) w[k][i]=wg[k]; else w[k][i]=0.0;}
2840 // end loop through alignment columns i
2841 ///////////////////////////////////////////////////////////////////////
2844 for (j=1; j<=L; j++){
2845 delete[](n[j]); (n[j]) = NULL;
2847 delete[](n); (n) = NULL;
2856 * take sequence data from Clustal and transfer it into
2857 * hhalign accessible information (structure/class)
2859 * Note that hhalign does not see all sequences/profiles
2860 * but only sequences that are elements of the 2 profiles
2863 * References to the required sequences are passed into hhalign
2864 * through auxilliary pointers that are shallow copies of the
2865 * sequence/profile data available to Clustal.
2867 * Re-allocating memory for these auxilliary pointers
2868 * would be desaterous, as it might detach the memory
2872 Alignment::Transfer(char **ppcProf, int iCnt){
2874 /* @<variables local to Transfer@> */
2875 int iLen; /* length of profile */
2876 int k; /* generic iterator */
2878 /* @<initialisation@> */
2880 N_filtered = N_ss = 0;
2881 kss_dssp = ksa_dssp = kss_pred = kss_conf = -1;
2883 strcpy(longname, "unknown_long_seq_name");
2884 strcpy(name, "unknown_seq_name");
2885 strcpy(file, "unknown_file_name");
2888 /* @<determine length of profile@>
2889 all sequences in profile should have same length,
2890 so only do it for 1st */
2891 for (iLen = 0; '\0' != ppcProf[0][iLen]; iLen++);
2893 /* @<allocate memory for sequences etc@> */
2894 for (k = 0; k < iCnt; k++){
2895 #define GOOD_MEASURE 1000 /* Temporary -- can be removed once rest in place */
2896 I[k] = new(short unsigned int[iLen+2+GOOD_MEASURE]);
2897 X[k] = new(char[iLen+2+GOOD_MEASURE]);
2898 seq[k] = new(char[iLen+2+GOOD_MEASURE]);
2901 if (NULL == ppcProf[k]){
2902 printf("%s:%d: Arena[%d]=NULL, cnt=%d\n", __FILE__, __LINE__, k, iCnt);
2905 strcat(seq[k], ppcProf[k]);
2906 keep[k] = KEEP_CONDITIONALLY;
2907 display[k] = KEEP_CONDITIONALLY;
2908 sname[k] = new(char[GOOD_MEASURE]);
2909 strcpy(sname[k], "unknown_sname");
2910 } /* (0 <= k < iCnt) */
2911 /* FIXME: Soeding always makes 1st sequence permanent */
2912 /*keep[0] = KEEP_ALWAYS;
2913 display[k] = KEEP_ALWAYS;*/
2915 /* Believe that the first and last positions are
2916 most important in stability of this algorithm.
2917 Must make sure that at least 2 sequences with
2918 residues in these positions are kept.
2919 Think any sequence will do, but better to keep
2920 the one with the longest 'contig'
2922 int iSeq; /* sequence iterator */
2923 int iHeadLen = 0, iHeadID = -1; /* length & ID of longest head contig */
2924 int iTailLen = 0, iTailID = -1; /* length & ID of longest head contig */
2926 char *pcFind = NULL;
2929 printf("%s:%s:%d: NEW PROFILE (%d seq) ================\n",
2930 __FUNCTION__, __FILE__, __LINE__, iCnt);
2932 for (iSeq = 0; iSeq < iCnt; iSeq++){
2934 printf("%s:%s:%d: consider seq %d ------------------\n",
2935 __FUNCTION__, __FILE__, __LINE__, iSeq);
2937 pcFind = strchr(&seq[iSeq][1], '-');
2938 if (NULL == pcFind){
2939 /* no gap at all in this sequences, spans entire profile */
2940 iHeadID = iTailID = iSeq;
2941 iHeadLen = iTailLen = iLen;
2944 iCont = (int)(pcFind - &seq[iSeq][1]);
2945 if (iCont > iHeadLen){
2949 pcFind = strrchr(seq[iSeq], '-');
2950 iCont = iLen - (int)(pcFind - seq[iSeq]);
2951 if (iCont > iTailLen){
2957 printf("%s:%s:%d: seq %3d: len = %d(%d) %s\n",
2958 __FUNCTION__, __FILE__, __LINE__, iSeq, iCont, iLen, seq[iSeq]);
2960 } /* 0 <= iSeq < iCnt */
2962 printf("%s:%s:%d: seq %d is winner with head contig of %d, seq %d tail contig of %d\n"
2963 , __FUNCTION__, __FILE__, __LINE__, iHeadID, iHeadLen, iTailID, iTailLen);
2965 if ( (-1 == iHeadID) || (-1 == iTailID) ){
2966 printf("%s:%s:%d: profile has no leading and/or trailing residues (h=%d:t=%d:#=%d)\n",
2967 __FUNCTION__, __FILE__, __LINE__, iHeadID, iTailID, iCnt);
2970 keep[iHeadID] = KEEP_ALWAYS;
2971 keep[iTailID] = KEEP_ALWAYS;
2977 } /* this is the end of Transfer() */
2981 /* @* Alignment::ClobberGlobal (eg: qali)
2983 * Note: originally hhalign() was stand-alone code,
2984 * there are a couple of GLOBAL (!) variables,
2985 * which would have been destroyed on exit.
2986 * However, now there is no 'exit' from hhalign(),
2987 * and on re-entry the global variable must be clean again.
2990 Alignment::ClobberGlobal(void){
2993 these are essential to re-set (as some of them are used as flags) */
2994 for(int k=0; k<N_in; k++)
2996 delete[] sname[k]; sname[k] = NULL;
2997 delete[] seq[k]; seq[k] = NULL;
2998 delete[] X[k]; X[k] = NULL;
2999 delete[] I[k]; I[k] = NULL;
3001 delete[] nres; nres = NULL;
3002 delete[] first; first = NULL;
3003 delete[] last; last = NULL;
3004 delete[] ksort; ksort = NULL;
3005 N_in = N_filtered = n_display = 0;
3007 kss_dssp = ksa_dssp = kss_pred = kss_conf = kfirst = -1;
3009 /* @<re-set but keep memory@>
3010 do not free the memory but re-set content */
3011 longname[0] = '\0'; //delete[] longname; longname = NULL;
3012 keep[0] = '\0'; //delete[] keep; keep = NULL;
3013 display[0] = '\0'; //delete[] display; display = NULL;
3014 wg[0] = 0; //delete[] wg; wg = NULL;
3015 nseqs[0] = 0; //delete[] nseqs; nseqs = NULL;
3019 //delete[] sname; sname = NULL;
3020 //delete[] seq; seq = NULL;
3021 //delete[] X; X = NULL;
3022 //delete[] I; I = NULL;
3023 //delete[] l; l = NULL;