binaries/src/clustalo/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 **/