← Back to branch summary

~ubuntu-branches/ubuntu/quantal/clustalo/quantal

« back to all changes in this revision

Viewing changes to src/hhalign/hhhalfalignment-C.h

  • Committer: Package Import Robot
  • Author(s): Olivier Sallou
  • Date: 2011-12-14 11:21:56 UTC
  • Revision ID: package-import@ubuntu.com-20111214112156-y3chwm0t4yn2nevm
Tags: upstream-1.0.3
Import upstream version 1.0.3

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/hhalign/hhhalfalignment-C.h
 
1
/* -*- mode: c; tab-width: 4; c-basic-offset: 4; indent-tabs-mode: nil -*- */
 
2
 
 
3
/*********************************************************************
 
4
 * Clustal Omega - Multiple sequence alignment
 
5
 *
 
6
 * Copyright (C) 2010 University College Dublin
 
7
 *
 
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.
 
12
 *
 
13
 * This file is part of Clustal-Omega.
 
14
 *
 
15
 ********************************************************************/
 
16
 
 
17
/*
 
18
 *  RCS $Id: hhhalfalignment-C.h 227 2011-03-28 17:03:09Z fabian $
 
19
 */
 
20
 
 
21
// hhfullalignment.C
 
22
 
 
23
#ifndef MAIN
 
24
#define MAIN
 
25
#include <iostream>   // cin, cout, cerr
 
26
#include <fstream>    // ofstream, ifstream 
 
27
#include <stdio.h>    // printf
 
28
#include <stdlib.h>   // exit
 
29
#include <string>     // strcmp, strstr
 
30
#include <math.h>     // sqrt, pow
 
31
#include <limits.h>   // INT_MIN
 
32
#include <float.h>    // FLT_MIN
 
33
#include <time.h>     // clock
 
34
#include <ctype.h>    // islower, isdigit etc
 
35
using std::ios;
 
36
using std::ifstream;
 
37
using std::ofstream;
 
38
using std::cout;
 
39
using std::cerr;
 
40
using std::endl;
 
41
#include "util-C.h"     // imax, fmax, iround, iceil, ifloor, strint, strscn, strcut, substr, uprstr, uprchr, Basename etc.
 
42
#include "list.h"     // list data structure
 
43
#include "hash.h"     // hash data structure
 
44
#include "hhdecl-C.h"      // constants, class 
 
45
#include "hhutil-C.h"      // imax, fmax, iround, iceil, ifloor, strint, strscn, strcut, substr, uprstr, uprchr, Basename etc.
 
46
#include "hhhmm.h"       // class HMM
 
47
#include "hhalignment.h" // class Alignment
 
48
#include "hhhit.h"
 
49
#endif
 
50
 
 
51
/////////////////////////////////////////////////////////////////////////////
 
52
/////////////////////////////////////////////////////////////////////////////
 
53
// Methods of class HalfAlignment
 
54
/////////////////////////////////////////////////////////////////////////////
 
55
/////////////////////////////////////////////////////////////////////////////
 
56
  
 
57
 
 
58
 
 
59
/////////////////////////////////////////////////////////////////////////////
 
60
// Constructor
 
61
HalfAlignment::HalfAlignment(int maxseqdis)
 
62
{
 
63
  n=0; 
 
64
  sname=seq=NULL; 
 
65
  nss_dssp = nss_pred = nss_conf = nsa_dssp = ncons= -1;
 
66
  h = new(int[maxseqdis]);   //h[k] = next position of sequence k to be written
 
67
  s = new(char*[maxseqdis]);  //s[k][h] = character in column h, sequence k of output alignment
 
68
  l = new(int*[maxseqdis]);   //counts non-gap residues: l[k][i] = index of last residue AT OR BEFORE match state i in seq k
 
69
  m = new(int*[maxseqdis]);   //counts positions:        m[k][i] = position of match state i in string seq[k]  
 
70
}
 
71
 
 
72
/////////////////////////////////////////////////////////////////////////////////////
 
73
// Destructor
 
74
HalfAlignment::~HalfAlignment()
 
75
{
 
76
  Unset();
 
77
  delete[] h; h = NULL;
 
78
  delete[] s; s = NULL;
 
79
  delete[] l; l = NULL;
 
80
  delete[] m; m = NULL;
 
81
}
 
82
 
 
83
 
 
84
 
 
85
/////////////////////////////////////////////////////////////////////////////////////
 
86
/**
 
87
 * @brief Free memory in HalfAlignment arrays s[][], l[][], and m[][]
 
88
 */
 
89
void 
 
90
HalfAlignment::Unset()
 
91
{
 
92
   // Free memory for alignment characters and residue counts
 
93
  for (int k=0; k<n; k++) 
 
94
    {
 
95
      delete[] s[k]; s[k] = NULL;
 
96
      delete[] l[k]; l[k] = NULL;
 
97
      delete[] m[k]; m[k] = NULL;
 
98
    }
 
99
  n=0; 
 
100
  sname=seq=NULL; 
 
101
  nss_dssp = nss_pred = nss_conf = nsa_dssp = ncons= -1;
 
102
}
 
103
 
 
104
 
 
105
//////////////////////////////////////////////////////////////////////////////
 
106
/**
 
107
 * @brief Prepare a2m/a3m alignment: 
 
108
 * Calculate l[k][i] (residue indices) and m[k][i] (position in seq[k]) 
 
109
*/
 
110
void 
 
111
HalfAlignment::Set(char* name, char** seq_in, char** sname_in, int n_in, int L_in, int n1, int n2, int n3, int n4, int nc, int L_in2/*<--FS*/)
 
112
{
 
113
    int i;     /* counts match states in seq[k] */
 
114
    int ll;    /* counts residues LEFT from or at current position in seq[k] */
 
115
    int mm;    /* counts postions in string seq[k] */
 
116
    int k;     /* counts sequences */
 
117
    char c;
 
118
    char warned=0;
 
119
    
 
120
    nss_dssp=n1; nss_pred=n2; nss_conf=n3; nsa_dssp=n4; ncons=nc;
 
121
    seq=seq_in;     /* flat copy of sequences */
 
122
    sname=sname_in; /* flat copy of sequence names */
 
123
    n=n_in;
 
124
    L=L_in;    
 
125
    pos=0;
 
126
 
 
127
    /* Allocate memory for alignment characters and residue counts */
 
128
    for (k=0; k<n; k++)  {
 
129
        s[k]=new char[LINELEN];
 
130
        l[k]=new int[L+10+L_in2/*<--FS*/]; 
 
131
        m[k]=new int[L+10+L_in2/*<--FS*/];
 
132
        if (!s[k] || !l[k] || !m[k]) MemoryError("space for formatting HMM-HMM alignment");
 
133
        h[k]=0; //starting positions in alignment = 0
 
134
    } /* k <= 0 < n (= n_in) */
 
135
 
 
136
  for (k=0; k<n; k++) {
 
137
      m[k][0]=0;  // 0'th match state (virtual) is begin state at 0
 
138
      //if k is consensus sequence 
 
139
      if (k==nc) {
 
140
          for (i=1; i<=L; i++) m[k][i]=l[k][i]=i; 
 
141
          m[k][L+1]=l[k][L+1]=L; 
 
142
          continue;
 
143
      }
 
144
      i=1; mm=1; ll=1;
 
145
      while ((c=seq[k][mm]))
 
146
          {
 
147
              if (MatchChr(c)==c)    //count match/delete states
 
148
                  {
 
149
                      l[k][i]=ll;
 
150
                      m[k][i]=mm;
 
151
                      i++;
 
152
                  }
 
153
              if (WordChr(c)) ll++;  //index of next residue
 
154
              mm++;
 
155
          }
 
156
      l[k][i]=ll-1; //set l[k][L+1] eq number of residues in seq k (-1 since there is no residue at L+1st match state)
 
157
      m[k][i]=mm;   //set m[k][L+1]
 
158
 
 
159
      if ((i-1)!=L && !warned) 
 
160
          {
 
161
              cerr<<"Warning: sequence "<<sname[k]<<" in HMM "<<name<<" has "<<i<<" match states but should have "<<L<<"\n"; 
 
162
              warned=1;
 
163
          }
 
164
  } /* k <= 0 < n (= n_in) */
 
165
  //DEBUG
 
166
  if (v>=5)
 
167
      {
 
168
          printf("  i chr   m   l\n");
 
169
          for(i=0;i<=L+1;i++) printf("%3i   %1c %3i %3i\n",i,seq[0][m[0][i]],m[0][i],l[0][i]);
 
170
          printf("\n");
 
171
      }
 
172
 
 
173
} /*** end HalfAlignment::Set() ***/
 
174
 
 
175
 
 
176
/////////////////////////////////////////////////////////////////////////////////////
 
177
/**
 
178
 * @brief Fill in insert states following match state i (without inserting '.' to fill up)
 
179
 */
 
180
void 
 
181
HalfAlignment::AddInserts(int i)
 
182
{
 
183
  for (int k=0; k<n; k++)                        // for all sequences...
 
184
    for (int mm=m[k][i]+1; mm<m[k][i+1]; mm++)   // for all inserts between match state i and i+1...
 
185
      s[k][h[k]++]=seq[k][mm];                   // fill inserts into output alignment s[k]
 
186
}
 
187
 
 
188
/////////////////////////////////////////////////////////////////////////////////////
 
189
/**
 
190
 * @brief Fill up alignment with gaps '.' to generate flush end (all h[k] equal)
 
191
 */
 
192
void 
 
193
HalfAlignment::FillUpGaps()
 
194
{ 
 
195
  int k;      //counts sequences
 
196
  pos=0;
 
197
 
 
198
  // Determine max position h[k]
 
199
  for (k=0; k<n; k++) pos = imax(h[k],pos);
 
200
  
 
201
  // Fill in gaps up to pos
 
202
  for (k=0; k<n; k++) 
 
203
    {
 
204
      for (int hh=h[k]; hh<pos; hh++) s[k][hh]='.';
 
205
      h[k]=pos;
 
206
    }
 
207
}
 
208
 
 
209
/////////////////////////////////////////////////////////////////////////////////////
 
210
/**
 
211
 * @brief Fill in insert states following match state i and fill up gaps with '.' 
 
212
 */
 
213
void 
 
214
HalfAlignment::AddInsertsAndFillUpGaps(int i) 
 
215
{ 
 
216
  AddInserts(i); 
 
217
  FillUpGaps(); 
 
218
}
 
219
 
 
220
/////////////////////////////////////////////////////////////////////////////////////
 
221
/**
 
222
 * @brief Add gap column '.'
 
223
 */
 
224
void 
 
225
HalfAlignment::AddChar(char c)  
 
226
{ 
 
227
  for (int k=0; k<n; k++) s[k][h[k]++]=c;                
 
228
  pos++;
 
229
}
 
230
 
 
231
/////////////////////////////////////////////////////////////////////////////////////
 
232
/**
 
233
 * @brief Add match state column i as is
 
234
 */
 
235
void 
 
236
HalfAlignment::AddColumn(int i) 
 
237
{ 
 
238
  for (int k=0; k<n; k++) s[k][h[k]++]=seq[k][m[k][i]];  
 
239
  pos++;
 
240
}
 
241
 
 
242
/////////////////////////////////////////////////////////////////////////////////////
 
243
/**
 
244
 * @brief Add match state column i as insert state
 
245
 */
 
246
void 
 
247
HalfAlignment::AddColumnAsInsert(int i) 
 
248
{ 
 
249
  char c; 
 
250
  for (int k=0; k<n; k++) 
 
251
    if ((c=seq[k][m[k][i]])!='-' && (c<'0' || c>'9')) 
 
252
      s[k][h[k]++]=InsertChr(c); 
 
253
  pos++;
 
254
}
 
255
 
 
256
/////////////////////////////////////////////////////////////////////////////////////
 
257
/**
 
258
 * @brief Remove all characters c from template sequences
 
259
 */
 
260
void 
 
261
HalfAlignment::RemoveChars(char c)
 
262
{ 
 
263
  int k,h,hh;
 
264
  for (k=0; k<n; k++)
 
265
    {
 
266
      for (h=hh=0; h<pos; h++)
 
267
        if (s[k][h]!=c) s[k][hh++]=s[k][h];
 
268
      s[k][++hh]='\0';
 
269
    }
 
270
}
 
271
 
 
272
/////////////////////////////////////////////////////////////////////////////////////
 
273
/**
 
274
 * @brief Transform alignment sequences from A3M to A2M (insert ".")
 
275
 */
 
276
void 
 
277
HalfAlignment::BuildFASTA()
 
278
{
 
279
  AddInserts(0); 
 
280
  FillUpGaps();
 
281
  for (int i=1; i<=L; i++)
 
282
    {
 
283
      AddColumn(i); 
 
284
      AddInserts(i); 
 
285
      FillUpGaps();
 
286
    }
 
287
  ToFASTA();
 
288
}
 
289
 
 
290
/////////////////////////////////////////////////////////////////////////////////////
 
291
/**
 
292
 * @brief Transform alignment sequences from A3M to A2M (insert ".")
 
293
 */
 
294
void 
 
295
HalfAlignment::BuildA2M()
 
296
{
 
297
  AddInserts(0); 
 
298
  FillUpGaps();
 
299
  for (int i=1; i<=L; i++)
 
300
    {
 
301
      AddColumn(i); 
 
302
      AddInserts(i); 
 
303
      FillUpGaps();
 
304
    }
 
305
  AddChar('\0');
 
306
}
 
307
 
 
308
/////////////////////////////////////////////////////////////////////////////////////
 
309
/**
 
310
 * @brief Transform alignment sequences from A3M to A2M (insert ".")
 
311
 */
 
312
void 
 
313
HalfAlignment::BuildA3M()
 
314
{
 
315
  AddInserts(0); 
 
316
  for (int i=1; i<=L; i++)
 
317
    {
 
318
      AddColumn(i); 
 
319
      AddInserts(i); 
 
320
    }
 
321
  AddChar('\0');
 
322
}
 
323
 
 
324
/////////////////////////////////////////////////////////////////////////////////////
 
325
/**
 
326
 * @brief Transform alignment sequences from A2M to FASTA ( lowercase to uppercase and '.' to '-')
 
327
 */
 
328
void 
 
329
HalfAlignment::ToFASTA()
 
330
{
 
331
  for (int k=0; k<n; k++)
 
332
    {
 
333
      uprstr(s[k]);
 
334
      strtr(s[k],".","-");
 
335
    }
 
336
}
 
337
 
 
338
/////////////////////////////////////////////////////////////////////////////////////
 
339
/**
 
340
 * @brief Align query (HalfAlignment) to template (i.e. hit) match state structure
 
341
 */
 
342
void 
 
343
HalfAlignment::AlignToTemplate(Hit& hit)
 
344
{
 
345
  int i,j;
 
346
  int step;    // column of the HMM-HMM alignment (first:nstep, last:1)
 
347
  char state;
 
348
 
 
349
  if(0) {  //par.loc==0) { //////////////////////////////////////////// STILL NEEDED??
 
350
    // If in global mode: Add part of alignment before first MM state
 
351
    AddInserts(0); // Fill in insert states before first match state
 
352
    for (i=1; i<hit.i[hit.nsteps]; i++)
 
353
      {
 
354
        AddColumnAsInsert(i);
 
355
        AddInserts(i);
 
356
        if (par.outformat<=2) FillUpGaps();
 
357
      }
 
358
  }
 
359
 
 
360
  // Add endgaps (First state must be an MM state!!)
 
361
  for (j=1; j<hit.j[hit.nsteps]; j++)    
 
362
    {
 
363
      AddChar('-');
 
364
    }
 
365
 
 
366
  // Add alignment between first and last MM state
 
367
  for (step=hit.nsteps; step>=1; step--) 
 
368
  {
 
369
    state = hit.states[step];
 
370
    i = hit.i[step];
 
371
 
 
372
    switch(state)
 
373
      {
 
374
      case MM:  //MM pair state (both query and template in Match state)
 
375
        AddColumn(i);
 
376
        AddInserts(i);
 
377
        break;
 
378
      case DG: //D- state
 
379
      case MI: //MI state
 
380
        AddColumnAsInsert(i);
 
381
        AddInserts(i);
 
382
        break;
 
383
      case GD: //-D state
 
384
      case IM: //IM state
 
385
        AddChar('-');
 
386
        break;
 
387
      }
 
388
    if (par.outformat<=2) FillUpGaps();
 
389
 
 
390
  }
 
391
 
 
392
  if(0) { //par.loc==0) { //////////////////////////////////////////// STILL NEEDED??
 
393
 
 
394
    // If in global mode: Add part of alignment after last MM state
 
395
    for (i=hit.i[1]+1; i<=L; i++)    
 
396
      {
 
397
        AddColumnAsInsert(i);
 
398
        AddInserts(i);
 
399
        if (par.outformat==2) FillUpGaps();
 
400
      }
 
401
  }
 
402
 
 
403
  // Add endgaps 
 
404
  for (j=hit.j[1]+1; j<=hit.L; j++)    
 
405
    {
 
406
      AddChar('-');
 
407
    }
 
408
 
 
409
  // Add end-of-string character
 
410
  AddChar('\0');
 
411
}
 
412
 
 
413
 
 
414
/////////////////////////////////////////////////////////////////////////////////////
 
415
/**
 
416
 * @brief Write the a2m/a3m alignment into alnfile 
 
417
 */
 
418
void 
 
419
HalfAlignment::Print(char* alnfile)
 
420
{
 
421
  int k;      //counts sequences
 
422
  int omitted=0; // counts number of sequences with no residues in match states
 
423
  FILE *outf;
 
424
  if (strcmp(alnfile,"stdout"))
 
425
    {
 
426
      if (par.append) outf=fopen(alnfile,"a"); else outf=fopen(alnfile,"w");
 
427
      if (!outf) OpenFileError(alnfile);
 
428
    } 
 
429
  else
 
430
    outf = stdout;
 
431
  if (v>=3) cout<<"Writing alignment to "<<alnfile<<"\n";
 
432
 
 
433
  for (k=0; k<n; k++)
 
434
    {
 
435
      // Print sequence only if it contains at least one residue in a match state
 
436
      if (1) //strpbrk(s[k],"ABCDEFGHIKLMNPQRSTUVWXYZ1234567890")) 
 
437
        {
 
438
          fprintf(outf,">%s\n",sname[k]);
 
439
          fprintf(outf,"%s\n",s[k]);
 
440
        } else {
 
441
          omitted++;
 
442
          if (v>=3) printf("%-14.14s contains no residue in match state. Omitting sequence\n",sname[k]);
 
443
        }
 
444
    }
 
445
  if (v>=2 && omitted) printf("Omitted %i sequences in %s which contained no residue in match state\n",omitted,alnfile);
 
446
  fclose(outf);
 
447
}
 
448
 
 
449
 
 
450
/** EOF hhhalfalignment-C.h **/