binaries/src/clustalo/src/hhalign/hhalign.cpp

   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: hhalign.cpp 309 2016-06-13 14:10:11Z fabian $
  19  */
  20
  21 /* hhalign.C:
  22  * Align a multiple alignment to an alignment or HMM
  23  * Print out aligned input sequences in a3m format
  24  * Compile:              g++ hhalign.C -o hhalign -I/usr/include/ -L/usr/lib -lpng -lz -O3 -fno-strict-aliasing
  25  * Compile with efence:  g++ hhalign.C -o hhalign -I/usr/include/ -lefence -L/usr/lib -lpng -lz -O -g
  26  *
  27  * Error codes: 0: ok  1: file format error  2: file access error
  28  *              3: memory error  4: internal numeric error
  29  *              5: command line error
  30  */
  31 #undef PNG           /* include options for making png files?
  32                         (will need the png library) */
  33 #define MAIN
  34 #include <iostream>   // cin, cout, cerr
  35 #include <fstream>    // ofstream, ifstream
  36 #include <string.h>     // strcmp, strstr
  37 #include <stdio.h>    // printf
  38 #include <stdlib.h>   // exit
  39 #include <math.h>     // sqrt, pow
  40 #include <limits.h>   // INT_MIN
  41 #include <float.h>    // FLT_MIN
  42 #include <time.h>     // clock
  43 #include <ctype.h>    // islower, isdigit etc
  44
  45 using std::cout;
  46 using std::cerr;
  47 using std::endl;
  48 using std::ios;
  49 using std::ifstream;
  50 using std::ofstream;
  51
  52 int iAux_GLOBAL;
  53
  54 #include "general.h"
  55 #include "util-C.h"     /* imax, fmax, iround, iceil, ifloor, strint, strscn,
  56                          strcut, substr, uprstr, uprchr, Basename etc. */
  57 #include "list-C.h"     // list data structure
  58 #include "hash-C.h"     // hash data structure
  59
  60 #include "hhdecl-C.h"      // Constants, global variables, struct Parameters
  61 #include "hhutil-C.h"      /* MatchChr, InsertChr, aa2i, i2aa, log2,
  62                             fast_log2, WriteToScreen, */
  63 #include "hhmatrices-C.h"  // BLOSUM50, GONNET, HSDM
  64 #include "hhhmm.h"       // class HMM
  65 #include "hhhit.h"       // class Hit
  66 #include "hhalignment.h" // class Alignment
  67 #include "hhhalfalignment.h" // class HalfAlignment
  68 #include "hhfullalignment.h" // class FullAlignment
  69 #include "hhhitlist.h"   // class Hit
  70 #include "hhhmm-C.h"       // class HMM
  71 #include "hhalignment-C.h" // class Alignment
  72 #include "hhhit-C.h"       // class Hit
  73 #include "hhhalfalignment-C.h" // class HalfAlignment
  74 #include "hhfullalignment-C.h" // class FullAlignment
  75 #include "hhhitlist-C.h"   // class HitList
  76 #include "hhfunc-C.h"      // some functions common to hh programs
  77
  78 #ifdef PNG
  79 #include "pngwriter.h"   //PNGWriter (http://pngwriter.sourceforge.net/)
  80 #include "pngwriter.cc"  //PNGWriter (http://pngwriter.sourceforge.net/)
  81 #endif
  82
  83 ////////////////////////////////////////////////////////////////////////
  84 // Global variables
  85 ////////////////////////////////////////////////////////////////////////
  86 HMM *q;                 // Create query HMM with maximum of MAXRES match states
  87 HMM *t;                 /* Create template HMM with maximum of MAXRES
  88                           match states  */
  89 Alignment *qali;        /* (query alignment might be needed outside of hhfunc.C
  90                           for -a option) */
  91 Hit hit;               // Ceate new hit object pointed at by hit
  92 HitList hitlist;       /* list of hits with one Hit object for each
  93                           pairwise comparison done */
  94 char aliindices[256];  /* hash containing indices of all alignments
  95                           which to show in dot plot */
  96 char* dmapfile=NULL;   /* where to write the coordinates for the HTML map file
  97                           (to click the alignments) */
  98 char* alitabfile=NULL; // where to write pairs of aligned residues
  99 char* strucfile=NULL;  // where to read structure scores
 100 char* pngfile=NULL;    // pointer to pngfile
 101 char* tcfile=NULL;     // TCoffee output file name
 102 float probmin_tc=0.05; /* 5% minimum posterior probability for printing
 103                           pairs of residues for TCoffee */
 104
 105 int dotW=10;           // average score of dot plot over window [i-W..i+W]
 106 float dotthr=0.5;      // score threshold for dot plot
 107 int dotscale=600;      // size scale of dotplot
 108 char dotali=0;         // show no alignments in dotplot
 109 float dotsat=0.3;      // saturation of grid and alignments in dot plot
 110 float pself=0.001;     // maximum p-value of 2nd and following self-alignments
 111 int Nstochali=0;       // # of stochastically traced alignments in dot plot
 112 float** Pstruc=NULL;   /* structure matrix which can be multiplied to prob
 113                           ratios from aa column comparisons in Forward() */
 114 float** Sstruc=NULL;   /* structure matrix which can be added to log odds
 115                           from aa column comparisons in Forward() */
 116 bool nucleomode = false;
 117
 118 ////////////////////////////////////////////////////////////////////////////
 119 // Help functions
 120 ////////////////////////////////////////////////////////////////////////////
 121 /**
 122  * @brief  general help function, not accessible in Clustal-Omega
 123  */
 124 #if 0
 125 void
 126 help()
 127 {
 128   printf("\n");
 129   printf("HHalign %s\n",VERSION_AND_DATE);
 130   printf("Align a query alignment/HMM to a template alignment/HMM by HMM-HMM alignment\n");
 131   printf("If only one alignment/HMM is given it is compared to itself and the best\n");
 132   printf("off-diagonal alignment plus all further non-overlapping alignments above \n");
 133   printf("significance threshold are shown.\n");
 134   printf("%s",REFERENCE);
 135   printf("%s",COPYRIGHT);
 136   printf("\n");
 137   printf("Usage: %s -i query [-t template] [options]  \n",program_name);
 138   printf(" -i <file>     input query alignment  (fasta/a2m/a3m) or HMM file (.hhm)\n");
 139   printf(" -t <file>     input template alignment (fasta/a2m/a3m) or HMM file (.hhm)\n");
 140 #ifdef PNG
 141   printf(" -png <file>   write dotplot into PNG-file (default=none)           \n");
 142 #endif
 143   printf("\n");
 144   printf("Output options:                                                           \n");
 145   printf(" -o <file>     write output alignment to file\n");
 146   printf(" -ofas <file>  write alignments in FASTA, A2M (-oa2m) or A3M (-oa3m) format   \n");
 147   printf(" -a <file>     write query alignment in a3m format to file (default=none)\n");
 148   printf(" -aa <file>    append query alignment in a3m format to file (default=none)\n");
 149   printf(" -atab <file>  write alignment as a table (with posteriors) to file (default=none)\n");
 150   printf(" -v <int>      verbose mode: 0:no screen output  1:only warings  2: verbose\n");
 151   printf(" -seq  [1,inf[ max. number of query/template sequences displayed  (def=%i)  \n",par.nseqdis);
 152   printf(" -nocons       don't show consensus sequence in alignments (default=show) \n");
 153   printf(" -nopred       don't show predicted 2ndary structure in alignments (default=show) \n");
 154   printf(" -nodssp       don't show DSSP 2ndary structure in alignments (default=show) \n");
 155   printf(" -aliw int     number of columns per line in alignment list (def=%i)\n",par.aliwidth);
 156   printf(" -P <float>    for self-comparison: max p-value of alignments (def=%.2g\n",pself);
 157   printf(" -p <float>    minimum probability in summary and alignment list (def=%G) \n",par.p);
 158   printf(" -E <float>    maximum E-value in summary and alignment list (def=%G)     \n",par.E);
 159   printf(" -Z <int>      maximum number of lines in summary hit list (def=%i)       \n",par.Z);
 160   printf(" -z <int>      minimum number of lines in summary hit list (def=%i)       \n",par.z);
 161   printf(" -B <int>      maximum number of alignments in alignment list (def=%i)    \n",par.B);
 162   printf(" -b <int>      minimum number of alignments in alignment list (def=%i)    \n",par.b);
 163   printf(" -rank int     specify rank of alignment to write with -a or -aa option (default=1)\n");
 164   printf("\n");
 165 #ifdef PNG
 166   printf("Dotplot options:\n");
 167   printf(" -dwin <int>   average score in dotplot over window [i-W..i+W] (def=%i)      \n",dotW);
 168   printf(" -dthr <float> score threshold for dotplot (default=%.2f)                    \n",dotthr);
 169   printf(" -dsca <int>   if value <= 20: size of dot plot unit box in pixels           \n");
 170   printf("               if value > 20: maximum dot plot size in pixels (default=%i)   \n",dotscale);
 171   printf(" -dali <list>  show alignments with indices in <list> in dot plot            \n");
 172   printf("               <list> = <index1> ... <indexN>  or  <list> = all              \n");
 173   printf("\n");
 174 #endif
 175   printf("Filter input alignment (options can be combined):                         \n");
 176   printf(" -id   [0,100] maximum pairwise sequence identity (%%) (def=%i)   \n",par.max_seqid);
 177   printf(" -diff [0,inf[ filter most diverse set of sequences, keeping at least this    \n");
 178   printf("               many sequences in each block of >50 columns (def=%i)\n",par.Ndiff);
 179   printf(" -cov  [0,100] minimum coverage with query (%%) (def=%i) \n",par.coverage);
 180   printf(" -qid  [0,100] minimum sequence identity with query (%%) (def=%i) \n",par.qid);
 181   printf(" -qsc  [0,100] minimum score per column with query  (def=%.1f)\n",par.qsc);
 182 //   printf(" -csc  [0,100] minimum score per column with core alignment (def=%-.2f)\n",par.coresc);
 183 //   printf(" -qscc [0,100] minimum score per column of core sequence with query (def=%-.2f)\n",par.qsc_core);
 184   printf("\n");
 185   printf("Input alignment format:                                                     \n");
 186   printf(" -M a2m        use A2M/A3M (default): upper case = Match; lower case = Insert;\n");
 187   printf("               '-' = Delete; '.' = gaps aligned to inserts (may be omitted)   \n");
 188   printf(" -M first      use FASTA: columns with residue in 1st sequence are match states\n");
 189   printf(" -M [0,100]    use FASTA: columns with fewer than X%% gaps are match states   \n");
 190   printf("\n");
 191   printf("HMM-HMM alignment options:                                                  \n");
 192   printf(" -glob/-loc    global or local alignment mode (def=local)         \n");
 193   printf(" -alt <int>    show up to this number of alternative alignments (def=%i)    \n",par.altali);
 194   printf(" -vit          use Viterbi algorithm for alignment instead of MAC algorithm \n");
 195   printf(" -mac          use Maximum Accuracy (MAC) alignment (default)  \n");
 196   printf(" -mact [0,1[   posterior probability threshold for MAC alignment (def=%.3f) \n",par.mact);
 197   printf("               A threshold value of 0.0 yields global alignments.\n");
 198   printf(" -sto <int>    use global stochastic sampling algorithm to sample this many alignments\n");
 199   printf(" -excl <range> exclude query positions from the alignment, e.g. '1-33,97-168'\n");
 200   printf(" -shift [-1,1] score offset (def=%-.3f)                                      \n",par.shift);
 201   printf(" -corr [0,1]   weight of term for pair correlations (def=%.2f)               \n",par.corr);
 202   printf(" -ssm  0-4     0:no ss scoring [default=%i]               \n",par.ssm);
 203   printf("               1:ss scoring after alignment                                  \n");
 204   printf("               2:ss scoring during alignment                                 \n");
 205   printf(" -ssw  [0,1]   weight of ss score  (def=%-.2f)                               \n",par.ssw);
 206   printf("\n");
 207   printf(" -def          read default options from ./.hhdefaults or <home>/.hhdefault. \n");
 208   printf("\n");
 209   printf("Example: %s -i T0187.a3m -t d1hz4a_.hhm -png T0187pdb.png \n",program_name);
 210   cout<<endl;
 211 //   printf("More help:                                                         \n");
 212 //   printf(" -h out        output options                                      \n");
 213 //   printf(" -h hmm        options for building HMM from multiple alignment    \n");
 214 //   printf(" -h gap        options for setting gap penalties                   \n");
 215 //   printf(" -h ali        options for HMM-HMM alignment                       \n");
 216 //   printf(" -h all        all options \n");
 217 }
 218 #endif
 219
 220 /**
 221  * @brief  helpt for output, not accessible in Clustal-Omega
 222  */
 223 #if 0
 224 void
 225 help_out()
 226 {
 227   printf("\n");
 228   printf("Output options: \n");
 229   printf(" -v            verbose mode (default: show only warnings)                 \n");
 230   printf(" -v 0          suppress all screen output                                 \n");
 231   printf(" -nocons       don't show consensus sequence in alignments (default=show) \n");
 232   printf(" -nopred       don't show predicted 2ndary structure in alignments (default=show) \n");
 233   printf(" -nodssp       don't show DSSP SS 2ndary structure in alignments (default=show) \n");
 234   printf(" -seq  [1,inf[ max. number of query/template sequences displayed  (def=%i)  \n",par.nseqdis);
 235   printf(" -aliw [40,..[ number of columns per line in alignment list (def=%i)\n",par.aliwidth);
 236   printf(" -P <float>    for self-comparison: max p-value of alignments (def=%.2g\n",pself);
 237   printf(" -p <float>    minimum probability in summary and alignment list (def=%G) \n",par.p);
 238   printf(" -E <float>    maximum E-value in summary and alignment list (def=%G)     \n",par.E);
 239   printf(" -Z <int>      maximum number of lines in summary hit list (def=%i)       \n",par.Z);
 240   printf(" -z <int>      minimum number of lines in summary hit list (def=%i)       \n",par.z);
 241   printf(" -B <int>      maximum number of alignments in alignment list (def=%i)    \n",par.B);
 242   printf(" -b <int>      minimum number of alignments in alignment list (def=%i)    \n",par.b);
 243   printf(" -rank int     specify rank of alignment to write with -a or -aa option (def=1)\n");
 244   printf(" -tc <file>    write a TCoffee library file for the pairwise comparison   \n");
 245   printf(" -tct [0,100]  min. probobability of residue pairs for TCoffee (def=%i%%)\n",iround(100*probmin_tc));
 246   printf("\n");
 247 #ifdef PNG
 248   printf("Dotplot options:\n");
 249   printf(" -dwin int     average score in dotplot over window [i-W..i+W] (def=%i)   \n",dotW);
 250   printf(" -dthr float   score threshold for dotplot (default=%.2f)                 \n",dotthr);
 251   printf(" -dsca int     size of dot plot box in pixels  (default=%i)               \n",dotscale);
 252   printf(" -dali <list>  show alignments with indices in <list> in dot plot\n");
 253   printf("               <list> = <index1> ... <indexN>  or  <list> = all              \n");
 254   printf(" -dmap <file>  print list of coordinates in png plot  \n");
 255 #endif
 256 }
 257 #endif
 258
 259 /**
 260  * @brief  help hit HMM options, not accessible in Clustal-Omega
 261  */
 262 #if 0
 263 void
 264 help_hmm()
 265 {
 266   printf("\n");
 267   printf("Options to filter input alignment (options can be combined):              \n");
 268   printf(" -id   [0,100] maximum pairwise sequence identity (%%) (def=%i)   \n",par.max_seqid);
 269   printf(" -diff [0,inf[ filter most diverse set of sequences, keeping at least this    \n");
 270   printf("               many sequences in each block of >50 columns (def=%i)\n",par.Ndiff);
 271   printf(" -cov  [0,100] minimum coverage with query (%%) (def=%i) \n",par.coverage);
 272   printf(" -qid  [0,100] minimum sequence identity with query (%%) (def=%i) \n",par.qid);
 273   printf(" -qsc  [0,100] minimum score per column with query  (def=%.1f)\n",par.qsc);
 274 //   printf(" -csc  [0,100] minimum score per column with core alignment (def=%-.2f)\n",par.coresc);
 275 //   printf(" -qscc [0,100] minimum score per column of core sequence with query (def=%-.2f)\n",par.qsc_core);
 276   printf("                                                                          \n");
 277   printf("HMM-building options:                                                     \n");
 278   printf(" -M a2m        use A2M/A3M (default): upper case = Match; lower case = Insert;\n");
 279   printf("               '-' = Delete; '.' = gaps aligned to inserts (may be omitted)   \n");
 280   printf(" -M first      use FASTA: columns with residue in 1st sequence are match states\n");
 281   printf(" -M [0,100]    use FASTA: columns with fewer than X%% gaps are match states   \n");
 282   printf(" -tags         do NOT neutralize His-, C-myc-, FLAG-tags, and \n");
 283   printf("               trypsin recognition sequence to background distribution    \n");
 284   printf("                                                                          \n");
 285   printf("Pseudocount options:                                                      \n");
 286   printf(" -Gonnet       use the Gonnet substitution matrix (default)               \n");
 287   printf(" -Blosum50     use the Blosum50 substitution matrix                       \n");
 288   printf(" -Blosum62     use the Blosum62 substitution matrix                       \n");
 289   printf(" -HSDM         use the structure-derived HSDM substitution matrix         \n");
 290   printf(" -pcm  0-3     Pseudocount mode (default=%-i)                             \n",par.pcm);
 291   printf("               tau = substitution matrix pseudocount admixture            \n");
 292   printf("               0: no pseudo counts:     tau = 0                           \n");
 293   printf("               1: constant              tau = a                           \n");
 294   printf("               2: divergence-dependent: tau = a/(1 + ((Neff-1)/b)^c)       \n");
 295   printf("                  Neff=( (Neff_q^d+Neff_t^d)/2 )^(1/d)                       \n");
 296   printf("                  Neff_q = av number of different AAs per column in query  \n");
 297   printf("               3: column-specific:      tau = \'2\' * (Neff(i)/Neff)^(w*Neff/20)\n");
 298   printf(" -pca  [0,1]   set a (overall admixture) (def=%-.1f)                      \n",par.pca);
 299   printf(" -pcb  [1,inf[ set b (threshold for Neff) (def=%-.1f)                      \n",par.pcb);
 300   printf(" -pcc  [0,3]   set c (extinction exponent for tau(Neff))  (def=%-.1f)      \n",par.pcc);
 301   printf(" -pcw  [0,3]   set w (weight of pos-specificity for pcs) (def=%-.1f)      \n",par.pcw);
 302 }
 303 #endif
 304
 305 /**
 306  * @brief  help with gap handling, not accessible in Clustal-Omega
 307  */
 308 #if 0
 309 void
 310 help_gap()
 311 {
 312   printf("\n");
 313   printf("Gap cost options:                                                         \n");
 314   printf(" -gapb [0,inf[ transition pseudocount admixture (def=%-.2f)               \n",par.gapb);
 315   printf(" -gapd [0,inf[ Transition pseudocount admixture for opening gap (default=%-.2f)\n",par.gapd);
 316   printf(" -gape [0,1.5] Transition pseudocount admixture for extending gap (def=%-.1f)\n",par.gape);
 317   printf(" -gapf ]0,inf] factor for increasing/reducing the gap open penalty for deletes (def=%-.2f)\n",par.gapf);
 318   printf(" -gapg ]0,inf] factor for increasing/reducing the gap open penalty for deletes (def=%-.2f)\n",par.gapg);
 319   printf(" -gaph ]0,inf] factor for increasing/reducing the gap extension penalty for deletes(def=%-.2f)\n",par.gaph);
 320   printf(" -gapi ]0,inf] factor for increasing/reducing the gap extension penalty for inserts(def=%-.2f)\n",par.gapi);
 321   printf(" -egq  [0,inf[ penalty (bits) for end gaps aligned to query residues (def=%-.2f)\n",par.egq);
 322   printf(" -egt  [0,inf[ penalty (bits) for end gaps aligned to template residues (def=%-.2f)\n",par.egt);
 323 }
 324 #endif
 325
 326 /**
 327  * @brief  help with alignment options, not accessible in Clustal-Omega
 328  */
 329 #if 0
 330 void
 331 help_ali()
 332 {
 333   printf("\n");
 334   printf("Alignment options:  \n");
 335   printf(" -glob/-loc    global or local alignment mode (def=global)                \n");
 336   printf(" -mac          use Maximum Accuracy (MAC) alignment instead of Viterbi\n");
 337   printf(" -mact [0,1]   posterior prob threshold for MAC alignment (def=%.3f)      \n",par.mact);
 338   printf(" -sto <int>    use global stochastic sampling algorithm to sample this many alignments\n");
 339   printf(" -sc   <int>   amino acid score         (tja: template HMM at column j) (def=%i)\n",par.columnscore);
 340   printf("        0      = log2 Sum(tja*qia/pa)   (pa: aa background frequencies)    \n");
 341   printf("        1      = log2 Sum(tja*qia/pqa)  (pqa = 1/2*(pa+ta) )               \n");
 342   printf("        2      = log2 Sum(tja*qia/ta)   (ta: av. aa freqs in template)     \n");
 343   printf("        3      = log2 Sum(tja*qia/qa)   (qa: av. aa freqs in query)        \n");
 344    printf(" -corr [0,1]   weight of term for pair correlations (def=%.2f)            \n",par.corr);
 345   printf(" -shift [-1,1] score offset (def=%-.3f)                                   \n",par.shift);
 346   printf(" -r            repeat identification: multiple hits not treated as independent\n");
 347   printf(" -ssm  0-2     0:no ss scoring [default=%i]               \n",par.ssm);
 348   printf("               1:ss scoring after alignment                               \n");
 349   printf("               2:ss scoring during alignment                              \n");
 350   printf(" -ssw  [0,1]   weight of ss score compared to column score (def=%-.2f)    \n",par.ssw);
 351   printf(" -ssa  [0,1]   ss confusion matrix = (1-ssa)*I + ssa*psipred-confusion-matrix [def=%-.2f)\n",par.ssa);
 352 }
 353 #endif
 354
 355 /**
 356  * @brief  general help menu, not accessible in Clustal-Omega
 357  */
 358 #if 0
 359 void
 360 help_all()
 361 {
 362   help();
 363   help_out();
 364   help_hmm();
 365   help_gap();
 366   help_ali();
 367   printf("\n");
 368   printf("Default options can be specified in './.hhdefaults' or '~/.hhdefaults'\n");
 369 }
 370 #endif
 371
 372 /////////////////////////////////////////////////////////////////////////////
 373 //// Processing input options from command line and .hhdefaults file
 374 /////////////////////////////////////////////////////////////////////////////
 375 /**
 376  * @brief  process arguments from commandline, not accessible from Clustal-Omega
 377  */
 378 #if 0
 379 void
 380 ProcessArguments(int argc, char** argv)
 381 {
 382   //Processing command line input
 383   for (int i=1; i<argc; i++)
 384     {
 385       if (v>=4) cout<<i<<"  "<<argv[i]<<endl; //PRINT
 386       if (!strcmp(argv[i],"-i"))
 387         {
 388           if (++i>=argc || argv[i][0]=='-')
 389             {help(); cerr<<endl<<"Error in "<<program_name<<": no query file following -i\n"; exit(4);}
 390           else strcpy(par.infile,argv[i]);
 391         }
 392       else if (!strcmp(argv[i],"-t"))
 393         {
 394           if (++i>=argc || argv[i][0]=='-')
 395             {help(); cerr<<endl<<"Error in "<<program_name<<": no template file following -d\n"; exit(4);}
 396           else strcpy(par.tfile,argv[i]); /* FS, ProcessArguments */
 397         }
 398       else if (!strcmp(argv[i],"-o"))
 399         {
 400           if (++i>=argc)
 401             {help(); cerr<<endl<<"Error in "<<program_name<<": no filename following -o\n"; exit(4);}
 402           else strcpy(par.outfile,argv[i]);
 403         }
 404       else if (!strcmp(argv[i],"-ofas"))
 405         {
 406           par.outformat=1;
 407           if (++i>=argc || argv[i][0]=='-')
 408             {help() ; cerr<<endl<<"Error in "<<program_name<<": no output file following -o\n"; exit(4);}
 409           else strcpy(par.pairwisealisfile,argv[i]);
 410         }
 411       else if (!strcmp(argv[i],"-oa2m"))
 412         {
 413           par.outformat=2;
 414           if (++i>=argc || argv[i][0]=='-')
 415             {help() ; cerr<<endl<<"Error in "<<program_name<<": no output file following -o\n"; exit(4);}
 416           else strcpy(par.pairwisealisfile,argv[i]);
 417         }
 418       else if (!strcmp(argv[i],"-oa3m"))
 419         {
 420           par.outformat=3;
 421           if (++i>=argc || argv[i][0]=='-')
 422             {help() ; cerr<<endl<<"Error in "<<program_name<<": no output file following -o\n"; exit(4);}
 423           else strcpy(par.pairwisealisfile,argv[i]);
 424         }
 425       else if (!strcmp(argv[i],"-rank") && (i<argc-1)) par.hitrank=atoi(argv[++i]);
 426       else if (!strcmp(argv[i],"-Oa3m"))
 427         {
 428           par.append=0;
 429           if (++i>=argc || argv[i][0]=='-')
 430             {help() ; cerr<<endl<<"Error in "<<program_name<<": no output file following -Oa3m\n"; exit(4);}
 431           else strcpy(par.alnfile,argv[i]);
 432         }
 433       else if (!strcmp(argv[i],"-Aa3m"))
 434         {
 435           par.append=1;
 436           if (++i>=argc || argv[i][0]=='-')
 437             {help() ; cerr<<endl<<"Error in "<<program_name<<": no output file following -Aa3m\n"; exit(4);}
 438           else strcpy(par.alnfile,argv[i]);
 439         }
 440       else if (!strcmp(argv[i],"-Ohhm"))
 441         {
 442           par.append=0;
 443           if (++i>=argc || argv[i][0]=='-')
 444             {help() ; cerr<<endl<<"Error in "<<program_name<<": no output file following -Ohhm\n"; exit(4);}
 445           else strcpy(par.hhmfile,argv[i]);
 446         }
 447       else if (!strcmp(argv[i],"-Ahhm"))
 448         {
 449           par.append=1;
 450           if (++i>=argc || argv[i][0]=='-')
 451             {help() ; cerr<<endl<<"Error in "<<program_name<<": no output file following -Ahhm\n"; exit(4);}
 452           else strcpy(par.hhmfile,argv[i]);
 453         }
 454       else if (!strcmp(argv[i],"-Opsi"))
 455         {
 456           par.append=0;
 457           if (++i>=argc || argv[i][0]=='-')
 458             {help() ; cerr<<endl<<"Error in "<<program_name<<": no output file following -Opsi\n"; exit(4);}
 459           else strcpy(par.psifile,argv[i]);
 460         }
 461       else if (!strcmp(argv[i],"-Apsi"))
 462         {
 463           par.append=1;
 464           if (++i>=argc || argv[i][0]=='-')
 465             {help() ; cerr<<endl<<"Error in "<<program_name<<": no output file following -Apsi\n"; exit(4);}
 466           else strcpy(par.psifile,argv[i]);
 467         }
 468       else if (!strcmp(argv[i],"-png"))
 469         {
 470           if (++i>=argc)
 471             {help(); cerr<<endl<<"Error in "<<program_name<<": no filename following -png\n"; exit(4);}
 472           else
 473             {
 474               pngfile = new(char[strlen(argv[i])+1]);
 475               strcpy(pngfile,argv[i]);
 476             }
 477         }
 478       else if (!strcmp(argv[i],"-Struc"))
 479         {
 480           if (++i>=argc || argv[i][0]=='-')
 481             {help(); cerr<<endl<<"Error in "<<program_name<<": no query file following -Struc\n"; exit(4);}
 482           else
 483             {
 484               strucfile = new(char[strlen(argv[i])+1]);
 485               strcpy(strucfile,argv[i]);
 486             }
 487         }
 488       else if (!strcmp(argv[i],"-atab") || !strcmp(argv[i],"-Aliout"))
 489         {
 490           if (++i>=argc || argv[i][0]=='-')
 491             {help(); cerr<<endl<<"Error in "<<program_name<<": no query file following -Struc\n"; exit(4);}
 492           else
 493             {
 494               alitabfile = new(char[strlen(argv[i])+1]);
 495               strcpy(alitabfile,argv[i]);
 496             }
 497         }
 498       else if (!strcmp(argv[i],"-tc"))
 499         {
 500           if (++i>=argc || argv[i][0]=='-')
 501             {help() ; cerr<<endl<<"Error in "<<program_name<<": no output file following -Opsi\n"; exit(4);}
 502           else
 503             {
 504               tcfile = new(char[strlen(argv[i])+1]);
 505               strcpy(tcfile,argv[i]);
 506             }
 507         }
 508       else if (!strcmp(argv[i],"-h")|| !strcmp(argv[i],"--help"))
 509         {
 510           if (++i>=argc) {help(); exit(0);}
 511           if (!strcmp(argv[i],"out")) {help_out(); exit(0);}
 512           if (!strcmp(argv[i],"hmm")) {help_hmm(); exit(0);}
 513           if (!strcmp(argv[i],"gap")) {help_gap(); exit(0);}
 514           if (!strcmp(argv[i],"ali")) {help_ali(); exit(0);}
 515           if (!strcmp(argv[i],"all")) {help_all(); exit(0);}
 516           else {help(); exit(0);}
 517         }
 518       else if (!strcmp(argv[i],"-excl"))
 519         {
 520           if (++i>=argc) {help(); exit(4);}
 521           par.exclstr = new(char[strlen(argv[i])+1]);
 522           strcpy(par.exclstr,argv[i]);
 523         }
 524       else if (!strcmp(argv[i],"-v") && (i<argc-1) && argv[i+1][0]!='-' ) v=atoi(argv[++i]);
 525       else if (!strcmp(argv[i],"-v"))  v=2;
 526       else if (!strcmp(argv[i],"-v0")) v=0;
 527       else if (!strcmp(argv[i],"-v1")) v=1;
 528       else if (!strcmp(argv[i],"-v2")) v=2;
 529       else if (!strcmp(argv[i],"-v3")) v=3;
 530       else if (!strcmp(argv[i],"-v4")) v=4;
 531       else if (!strcmp(argv[i],"-v5")) v=5;
 532       else if (!strcmp(argv[i],"-P") && (i<argc-1)) pself=atof(argv[++i]);
 533       else if (!strcmp(argv[i],"-p") && (i<argc-1)) par.p = atof(argv[++i]);
 534       else if (!strcmp(argv[i],"-e") && (i<argc-1)) par.E = atof(argv[++i]);
 535       else if (!strcmp(argv[i],"-E") && (i<argc-1)) par.E = atof(argv[++i]);
 536       else if (!strcmp(argv[i],"-b") && (i<argc-1)) par.b = atoi(argv[++i]);
 537       else if (!strcmp(argv[i],"-B") && (i<argc-1)) par.B = atoi(argv[++i]);
 538       else if (!strcmp(argv[i],"-z") && (i<argc-1)) par.z = atoi(argv[++i]);
 539       else if (!strcmp(argv[i],"-Z") && (i<argc-1)) par.Z = atoi(argv[++i]);
 540       else if (!strncmp(argv[i],"-nocons",7)) par.showcons=0;
 541       else if (!strncmp(argv[i],"-nopred",7)) par.showpred=0;
 542       else if (!strncmp(argv[i],"-nodssp",7)) par.showdssp=0;
 543       else if (!strncmp(argv[i],"-mark",7)) par.mark=1;
 544       else if (!strcmp(argv[i],"-seq") && (i<argc-1))  par.nseqdis=atoi(argv[++i]);
 545       else if (!strcmp(argv[i],"-aliw") && (i<argc-1)) par.aliwidth=atoi(argv[++i]);
 546       else if (!strcmp(argv[i],"-id") && (i<argc-1))   par.max_seqid=atoi(argv[++i]);
 547       else if (!strcmp(argv[i],"-tct") && (i<argc-1))  probmin_tc=atoi(argv[++i]);
 548       else if (!strcmp(argv[i],"-dwin") && (i<argc-1)) dotW=atoi(argv[++i]);
 549       else if (!strcmp(argv[i],"-dsca") && (i<argc-1)) dotscale=atoi(argv[++i]);
 550       else if (!strcmp(argv[i],"-dthr") && (i<argc-1)) dotthr=atof(argv[++i]);
 551       else if (!strcmp(argv[i],"-dali") && (i<argc-1))
 552         {
 553           dotali=1;
 554           for (int index=0; index<256; index++) aliindices[index]=0;
 555           while (i+1<argc && argv[i+1][0]!='-') // adds index to hash aliindices
 556             {
 557               i++;
 558               if (strcmp(argv[i],"all")) aliindices[atoi(argv[i])]=1;
 559               else dotali=2;
 560             }
 561         }
 562       else if (!strcmp(argv[i],"-dmap"))
 563         {
 564           if (++i>=argc)
 565             {help(); cerr<<endl<<"Error in "<<program_name<<": no filename following -o\n"; exit(4);}
 566           else
 567             {
 568               dmapfile = new(char[strlen(argv[i])+1]);
 569               strcpy(dmapfile,argv[i]);
 570             }
 571         }
 572       else if (!strcmp(argv[i],"-dsat") && (i<argc-1)) dotsat=atof(argv[++i]);
 573       else if (!strcmp(argv[i],"-qid") && (i<argc-1))  par.qid=atoi(argv[++i]);
 574       else if (!strcmp(argv[i],"-qsc") && (i<argc-1))  par.qsc=atof(argv[++i]);
 575       else if (!strcmp(argv[i],"-cov") && (i<argc-1))  par.coverage=atoi(argv[++i]);
 576       else if (!strcmp(argv[i],"-diff") && (i<argc-1)) par.Ndiff=atoi(argv[++i]);
 577       else if (!strcmp(argv[i],"-qscc") && (i<argc-1))    par.qsc_core=atof(argv[++i]);
 578       else if (!strcmp(argv[i],"-csc") && (i<argc-1))     par.coresc=atof(argv[++i]);
 579       else if (!strcmp(argv[i],"-Gonnet")) par.matrix=0;
 580       else if (!strcmp(argv[i],"-HSDM")) par.matrix=1;
 581       else if (!strcmp(argv[i],"-BLOSUM50")) par.matrix=2;
 582       else if (!strcmp(argv[i],"-Blosum50")) par.matrix=2;
 583       else if (!strcmp(argv[i],"-B50")) par.matrix=2;
 584       else if (!strcmp(argv[i],"-BLOSUM62")) par.matrix=3;
 585       else if (!strcmp(argv[i],"-Blosum62")) par.matrix=3;
 586       else if (!strcmp(argv[i],"-B62")) par.matrix=3;
 587       else if (!strcmp(argv[i],"-pcm") && (i<argc-1)) par.pcm=atoi(argv[++i]);
 588       else if (!strcmp(argv[i],"-pca") && (i<argc-1)) par.pca=atof(argv[++i]);
 589       else if (!strcmp(argv[i],"-pcb") && (i<argc-1)) par.pcb=atof(argv[++i]);
 590       else if (!strcmp(argv[i],"-pcc") && (i<argc-1)) par.pcc=atof(argv[++i]);
 591       else if (!strcmp(argv[i],"-pcw") && (i<argc-1)) par.pcw=atof(argv[++i]);
 592       else if (!strcmp(argv[i],"-gapb") && (i<argc-1)) { par.gapb=atof(argv[++i]); if (par.gapb<=0.01) par.gapb=0.01;}
 593       else if (!strcmp(argv[i],"-gapd") && (i<argc-1)) par.gapd=atof(argv[++i]);
 594       else if (!strcmp(argv[i],"-gape") && (i<argc-1)) par.gape=atof(argv[++i]);
 595       else if (!strcmp(argv[i],"-gapf") && (i<argc-1)) par.gapf=atof(argv[++i]);
 596       else if (!strcmp(argv[i],"-gapg") && (i<argc-1)) par.gapg=atof(argv[++i]);
 597       else if (!strcmp(argv[i],"-gaph") && (i<argc-1)) par.gaph=atof(argv[++i]);
 598       else if (!strcmp(argv[i],"-gapi") && (i<argc-1)) par.gapi=atof(argv[++i]);
 599       else if (!strcmp(argv[i],"-egq") && (i<argc-1)) par.egq=atof(argv[++i]);
 600       else if (!strcmp(argv[i],"-egt") && (i<argc-1)) par.egt=atof(argv[++i]);
 601       else if (!strcmp(argv[i],"-ssgap")) par.ssgap=1;
 602       else if (!strcmp(argv[i],"-ssgapd") && (i<argc-1)) par.ssgapd=atof(argv[++i]);
 603       else if (!strcmp(argv[i],"-ssgape") && (i<argc-1)) par.ssgape=atof(argv[++i]);
 604       else if (!strcmp(argv[i],"-ssgapi") && (i<argc-1)) par.ssgapi=atoi(argv[++i]);
 605       else if (!strcmp(argv[i],"-ssm") && (i<argc-1)) par.ssm=atoi(argv[++i]);
 606       else if (!strcmp(argv[i],"-ssw") && (i<argc-1)) par.ssw=atof(argv[++i]);
 607       else if (!strcmp(argv[i],"-ssa") && (i<argc-1)) par.ssa=atof(argv[++i]);
 608       else if (!strncmp(argv[i],"-glo",3)) {par.loc=0; if (par.mact>0.3 && par.mact<0.301) {par.mact=0;} }
 609       else if (!strncmp(argv[i],"-loc",3)) par.loc=1;
 610       else if (!strncmp(argv[i],"-alt",4) && (i<argc-1)) par.altali=atoi(argv[++i]);
 611       else if (!strcmp(argv[i],"-map") || !strcmp(argv[i],"-MAP")) par.forward=2;
 612       else if (!strcmp(argv[i],"-mac") || !strcmp(argv[i],"-MAC")) par.forward=2;
 613       else if (!strcmp(argv[i],"-vit")) par.forward=0;
 614       else if (!strcmp(argv[i],"-sto") && (i<argc-1))  {Nstochali=atoi(argv[++i]); par.forward=1;}
 615       else if (!strcmp(argv[i],"-r")) par.repmode=1;
 616       else if (!strcmp(argv[i],"-M") && (i<argc-1))
 617         if (!strcmp(argv[++i],"a2m") || !strcmp(argv[i],"a3m"))  par.M=1;
 618         else if(!strcmp(argv[i],"first"))  par.M=3;
 619         else if (argv[i][0]>='0' && argv[i][0]<='9') {par.Mgaps=atoi(argv[i]); par.M=2;}
 620         else cerr<<endl<<"WARNING: Ignoring unknown argument: -M "<<argv[i]<<"\n";
 621       else if (!strcmp(argv[i],"-shift") && (i<argc-1)) par.shift=atof(argv[++i]);
 622       else if (!strcmp(argv[i],"-mact") && (i<argc-1)) {par.mact=atof(argv[++i]); par.forward=2;}
 623       else if (!strcmp(argv[i],"-wstruc") && (i<argc-1)) par.wstruc=atof(argv[++i]);
 624       else if (!strcmp(argv[i],"-opt") && (i<argc-1)) par.opt=atoi(argv[++i]);
 625       else if (!strcmp(argv[i],"-sc") && (i<argc-1)) par.columnscore=atoi(argv[++i]);
 626       else if (!strcmp(argv[i],"-corr") && (i<argc-1)) par.corr=atof(argv[++i]);
 627       else if (!strcmp(argv[i],"-def")) par.readdefaultsfile=1;
 628       else if (!strcmp(argv[i],"-ovlp") && (i<argc-1)) par.min_overlap=atoi(argv[++i]);
 629       else if (!strcmp(argv[i],"-tags")) par.notags=0;
 630       else if (!strcmp(argv[i],"-notags")) par.notags=1;
 631       else cerr<<endl<<"WARNING: Ignoring unknown option "<<argv[i]<<" ...\n";
 632       if (v>=4) cout<<i<<"  "<<argv[i]<<endl; //PRINT
 633     } // end of for-loop for command line input
 634 }
 635 #endif
 636
 637
 638
 639
 640 /////////////////////////////////////////////////////////////////////////////
 641 //// MAIN PROGRAM
 642 /////////////////////////////////////////////////////////////////////////////
 643 /**
 644  *
 645  * hhalign()
 646  *
 647  * @param[in,out] ppcFirstProf
 648  * first profile to be aligned
 649  * @param[in] iFirstCnt
 650  * number of sequences in 1st profile
 651  * @param[in,out] ppcSecndProf
 652  * second profile to be aligned
 653  * @param[in] iSecndCnt
 654  * number of sequences in 2nd profile
 655  * @param[out] dScore_p
 656  * score of alignment
 657  * @param[in] prHMM
 658  * HMM info of external HMM (background)
 659  * @param[in] pcPrealigned1
 660  * (gapped) 1st sequence aligned to background HMM
 661  * @param[in] pcRepresent1
 662  * (gapped) sequence representing background HMM aligned to 1st sequence
 663  * @param[in] pcPrealigned2
 664  * (gapped) 2nd sequence aligned to background HMM
 665  * @param[in] pcRepresent2
 666  * (gapped) sequence representing background HMM aligned to 2nd sequence
 667  * @param[in] rHhalignPara,
 668  * various parameters passed down from commandline
 669  * e.g., iMaxRamMB
 670  * @param[out] rHHscores
 671  * qualify goodness of alignment
 672  *
 673  * iFlag,zcAux,zcError are debugging arguments
 674  *
 675  * @return Non-zero on error
 676  */
 677 extern "C" int
 678 hhalign(char **ppcFirstProf, int iFirstCnt, double *pdWeightsL,
 679         char **ppcSecndProf, int iSecndCnt, double *pdWeightsR,
 680         double *dScore_p, hmm_light *prHMM, hmm_light *prHMM2,
 681         char *pcPrealigned1, char *pcRepresent1,
 682         char *pcPrealigned2, char *pcRepresent2,
 683         hhalign_para rHhalignPara, hhalign_scores *rHHscores,
 684         int iFlag, int iVerbosity,
 685         char zcAux[], char zcError[]) {
 686
 687
 688 #ifdef CLUSTALO
 689     int iRetVal = RETURN_OK;
 690     iAux_GLOBAL = iFlag;
 691 #ifndef CLUSTALO_NOFILE
 692     int argc = 0;
 693     char **argv = NULL;
 694     argv = (char **)malloc(argc*sizeof(char *));
 695     for (int i = 0; i < argc; i++){
 696         argv[i] = (char *)malloc(100);
 697     }
 698     strcpy(argv[0], "./hhalign");
 699     strcpy(argv[1], "-t");
 700     strcpy(argv[2], "hhalign.C");
 701     strcpy(argv[3], "-i");
 702     strcpy(argv[4], "hhalign.C");
 703     strcpy(argv[5], "-ofas");
 704     strcpy(argv[6], "out");
 705 #endif
 706 #endif
 707     /*char* argv_conf[MAXOPT];*/ /* Input arguments from .hhdefaults file
 708       (first=1: argv_conf[0] is not used) */
 709     /*int argc_conf;*/           // Number of arguments in argv_conf
 710     /*char inext[IDLEN]="";*/    // Extension of query input file (hhm or a3m)
 711     /*char text[IDLEN]="";*/     // Extension of template input file (hhm or a3m)
 712     /*int** ali=NULL;*/          // ali[i][j]=1 if (i,j) is part of an alignment
 713     /*int** alisto=NULL;*/       // ali[i][j]=1 if (i,j) is part of an alignment
 714     int Nali;                    /* number of normally backtraced alignments
 715                                     in dot plot */
 716
 717     SetDefaults(rHhalignPara);
 718     strcpy(par.tfile,""); /* FS, Initialise Argument */
 719     strcpy(par.alnfile,"");
 720     par.p=0.0 ;                  /* minimum threshold for inclusion in hit list
 721                                     and alignment listing */
 722     par.E=1e6;                   /* maximum threshold for inclusion in hit list
 723                                     and alignment listing */
 724     par.b=1;                     // min number of alignments
 725     par.B=100;                   // max number of alignments
 726     par.z=1;                     // min number of lines in hit list
 727     par.Z=100;                   // max number of lines in hit list
 728     par.append=0;                /* append alignment to output file
 729                                     with -a option */
 730     par.altali=1;                /* find only ONE (possibly overlapping)
 731                                     subalignment  */
 732     par.hitrank=0;               /* rank of hit to be printed as a3m alignment
 733                                     (default=0) */
 734     par.outformat=3;             // default output format for alignment is a3m
 735     hit.self=0;                  // no self-alignment
 736     par.forward=2;               /* 0: Viterbi algorithm;
 737                                     1: Viterbi+stochastic sampling;
 738                                     2:Maximum Accuracy (MAC) algorithm */
 739
 740     // Make command line input globally available
 741 #ifndef CLUSTALO_NOFILE
 742     par.argv=argv;
 743     par.argc=argc;
 744     RemovePathAndExtension(program_name,argv[0]);
 745 #endif
 746
 747 #ifndef CLUSTALO_NOFILE
 748     /* Enable changing verbose mode before defaults file
 749        and command line are processed */
 750     for (int i=1; i<argc; i++)
 751         {
 752             if (!strcmp(argv[i],"-def")) par.readdefaultsfile=1;
 753             else if (argc>1 && !strcmp(argv[i],"-v0")) v=0;
 754             else if (argc>1 && !strcmp(argv[i],"-v1")) v=1;
 755             else if (argc>2 && !strcmp(argv[i],"-v")) v=atoi(argv[i+1]);
 756         }
 757
 758     // Read .hhdefaults file?
 759     if (par.readdefaultsfile)
 760         {
 761             // Process default otpions from .hhconfig file
 762             ReadDefaultsFile(argc_conf,argv_conf);
 763             ProcessArguments(argc_conf,argv_conf);
 764         }
 765 #endif
 766     /* Process command line options
 767        (they override defaults from .hhdefaults file) */
 768 #ifndef CLUSTALO_NOFILE
 769     ProcessArguments(argc,argv);
 770 #endif
 771
 772
 773 #ifdef CLUSTALO
 774     int iAuxLen1 = strlen(ppcFirstProf[0]);
 775     int iAuxLen2 = strlen(ppcSecndProf[0]);
 776     if ( (0 == iAuxLen1) || (0 == iAuxLen2) ){ /* problem with empty profiles, FS, r249 -> r250 */
 777         sprintf(zcError, "%s:%s:%d: strlen(prof1)=%d, strlen(prof2)=%d -- Nothing to align!\n",
 778                 __FUNCTION__, __FILE__, __LINE__, iAuxLen1, iAuxLen2);
 779         iRetVal = RETURN_UNKNOWN;
 780         /* Note: at this stage cannot do   'goto this_is_the_end;'
 781            because would cross initialisation of several variables   */
 782         return iRetVal;
 783     }
 784     par.maxResLen  = iAuxLen1 > iAuxLen2 ? iAuxLen1 : iAuxLen2;
 785     const int ciGoodMeasureSeq = 10;
 786     par.maxResLen += ciGoodMeasureSeq;
 787     par.maxColCnt  = iAuxLen1 + iAuxLen2 + ciGoodMeasureSeq;
 788     //par.max_seqid=iFirstCnt+iSecndCnt+3; /* -id     */
 789     par.max_seqid=DEFAULT_FILTER;        /* -id     */
 790     par.loc=0; par.mact=0;              /* -glob   */
 791     par.nseqdis=iFirstCnt+iSecndCnt;   /* -seq    */
 792     par.showcons=0;                   /* -nocons */
 793     par.showdssp=0;                  /* -nodssp */
 794     par.Mgaps=100;                  /* -M      */
 795     par.M=2;                       /* -M      */
 796     par.pdWg1=pdWeightsL;         /* tree wg */
 797     par.pdWg2=pdWeightsR;        /* tree wg */
 798     v = 0;                      /* -v0     */
 799     /* NOTE: *qali declared globally but only created here,
 800        pass in number of sequences to get rid of statically
 801        defined MAXSEQ (FS)
 802     */
 803     Alignment qali(iFirstCnt+iSecndCnt);
 804     HMM q(iFirstCnt+iSecndCnt);
 805     HMM t(iFirstCnt+iSecndCnt);
 806 #endif
 807
 808
 809 #ifndef CLUSTALO_NOFILE
 810     // Check command line input and default values
 811     if (!*par.infile)
 812         {help(); cerr<<endl<<"Error in "<<program_name<<": no query alignment file given (-i file)\n"; exit(4);}
 813     if (par.forward==2 && par.loc==0)
 814         {
 815             if (par.mact<0.301 || par.mact>0.300)
 816                 if (v>=1) fprintf(stderr,"REMARK: in -mac -global mode -mact is forced to 0\n");
 817             par.mact=0;
 818             //       par.loc=1; // global forward-backward algorithm seems to work fine! use it in place of local version.
 819         }
 820
 821     // Get rootname (no directory path, no extension) and extension of infile
 822     RemoveExtension(q.file,par.infile);
 823     RemoveExtension(t.file,par.tfile); /* FS, NOFILE2 (commented out) */
 824     Extension(inext,par.infile);
 825     Extension(text,par.tfile); /*FS, NOFILE2 (commented out) */
 826
 827     // Check option compatibilities
 828     if (par.nseqdis>MAXSEQDIS-3-par.showcons) par.nseqdis=MAXSEQDIS-3-par.showcons; //3 reserved for secondary structure
 829     if (par.aliwidth<20) par.aliwidth=20;
 830     if (par.pca<0.001) par.pca=0.001; // to avoid log(0)
 831     if (par.b>par.B) par.B=par.b;
 832     if (par.z>par.Z) par.Z=par.z;
 833     if (par.hitrank>0) par.altali=0;
 834
 835     // Input parameters
 836     if (v>=3)
 837         {
 838             cout<<"query file : "<<par.infile<<"\n";
 839             cout<<"template file: "<<par.tfile<<"\n";/*FS, NOFILE2 (commented out) */
 840             cout<<"Output file:  "<<par.outfile<<"\n";
 841             cout<<"Alignment file:  "<<par.alnfile<<"\n";
 842         }
 843 #endif /* NOFILE2 */
 844
 845
 846     // Set (global variable) substitution matrix and derived matrices
 847         // DD: new experimental matrices and params for nucleotides
 848         if(rHhalignPara.bIsDna)
 849         {
 850                 nucleomode = true;
 851                 SetDnaDefaults(rHhalignPara);
 852                 SetDnaSubstitutionMatrix();
 853         }
 854         else if(rHhalignPara.bIsRna)
 855         {
 856                 nucleomode = true;
 857                 SetRnaDefaults(rHhalignPara);
 858                 SetRnaSubstitutionMatrix();
 859         }
 860         else
 861                 SetSubstitutionMatrix();
 862
 863     // Set secondary structure substitution matrix
 864     SetSecStrucSubstitutionMatrix();
 865
 866
 867     /* moved Viterbi switch from after RnP() to here,
 868        switch after RnP() ineffectual as RnP decides log/lin of transition,
 869        however log/lin of transitions depends on MAC/Viterbi,
 870        FS, r228 -> r229 */
 871     int qL, tL;
 872     if (iFirstCnt > 0) {
 873         qL = strlen(ppcFirstProf[0]);
 874     }
 875     else {
 876         qL = prHMM->L;
 877     }
 878     if (iSecndCnt > 0) {
 879         tL = strlen(ppcSecndProf[0]);
 880     }
 881     else {
 882         tL = prHMM2->L;
 883     }
 884
 885
 886     //  const float MEMSPACE_DYNPROG = 512*1024*1024;
 887     /* determine amount of memory available for MAC on command-line; FS, r240 -> r241 */
 888     const float MEMSPACE_DYNPROG = (double)1024*1024*rHhalignPara.iMacRamMB;
 889     // longest allowable length of database HMM
 890     int Lmaxmem=(int)((float)MEMSPACE_DYNPROG/qL/6/8);
 891     if (par.forward==2 && tL+2>=Lmaxmem) {
 892         sprintf(zcError, "%s:%s:%d: switch to Viterbi (qL=%d, tL=%d, MAC-RAM=%d\n)",
 893                 __FUNCTION__, __FILE__, __LINE__, qL, tL, rHhalignPara.iMacRamMB);
 894         if (v>=1)
 895             cerr<<"WARNING: Not sufficient memory to realign with MAC algorithm. Using Viterbi algorithm."<<endl;
 896         par.forward=0;
 897
 898         /* use different 'fudge' parameters for Viterbi, FS, r228 -> r229 */
 899         par.pca = par.pcaV;
 900         par.pcb = par.pcbV;
 901         par.pcc = par.pccV;
 902         par.pcw = par.pcwV;
 903
 904         par.gapb = par.gapbV;
 905         par.gapd = par.gapdV;
 906         par.gape = par.gapeV;
 907         par.gapf = par.gapfV;
 908         par.gapg = par.gapgV;
 909         par.gaph = par.gaphV;
 910         par.gapi = par.gapiV;
 911
 912     }
 913
 914     // Read input file (HMM, HHM, or alignment format), and add pseudocounts etc.
 915     q.cQT = 'q';
 916     if (OK != ReadAndPrepare(INTERN_ALN_2_HMM,
 917                              ppcFirstProf, iFirstCnt, prHMM,
 918                              pcPrealigned1, pcRepresent1, pdWeightsL,
 919                              (char*)(""), q, &qali)) {
 920         sprintf(zcError, "%s:%s:%d: Problem Reading/Preparing q-profile (len=%d)\n",
 921                 __FUNCTION__, __FILE__, __LINE__, qL);
 922         iRetVal = RETURN_FROM_RNP;
 923         goto this_is_the_end;
 924     }
 925     // Set query columns in His-tags etc to Null model distribution
 926     if (par.notags) q.NeutralizeTags();
 927
 928     // Do self-comparison?
 929     if (0 /* !*par.tfile // FS, 2010-03-10 */)
 930         {
 931             // Deep-copy q into t
 932             t = q;
 933
 934             // Find overlapping alternative alignments
 935             hit.self=1;
 936         }
 937     // Read template alignment/HMM t and add pseudocounts
 938     else
 939         {
 940             char infile[] = "";
 941             /* Read input file (HMM, HHM, or alignment format),
 942                and add pseudocounts etc. */
 943             t.cQT = 't';
 944             if (OK != ReadAndPrepare(INTERN_ALN_2_HMM,
 945                                      ppcSecndProf, iSecndCnt, prHMM2,
 946                                      pcPrealigned2, pcRepresent2, pdWeightsR,
 947                                      infile, t)) {
 948                 sprintf(zcError, "%s:%s:%d: Problem Reading/Preparing t-profile (len=%d)\n",
 949                         __FUNCTION__, __FILE__, __LINE__, tL);
 950                 iRetVal = RETURN_FROM_RNP;
 951                 goto this_is_the_end;
 952             }
 953         }
 954
 955   // Factor Null model into HMM t
 956   t.IncludeNullModelInHMM(q,t);
 957
 958   /* alignment will fail if one profile does not contain useful characters, FS, r259 -> r260 */
 959   if ( (q.L <= 0) || (t.L <= 0) ){
 960       sprintf(zcError, "%s:%s:%d: Problem Reading/Preparing profiles (len(q)=%d/len(t)=%d)\n",
 961               __FUNCTION__, __FILE__, __LINE__, q.L, t.L);
 962       iRetVal = RETURN_FROM_RNP;
 963       goto this_is_the_end;
 964   }
 965
 966   /* switch at this stage is ineffectual as log/lin already decided in RnP().
 967      FS, r228 -> r229 */
 968   /*const float MEMSPACE_DYNPROG = 512*1024*1024;
 969   // longest allowable length of database HMM
 970   int Lmaxmem=(int)((float)MEMSPACE_DYNPROG/q.L/6/8);
 971   if (par.forward==2 && t.L+2>=Lmaxmem)
 972   {
 973   if (v>=1)
 974   cerr<<"WARNING: Not sufficient memory to realign with MAC algorithm. Using Viterbi algorithm."<<endl;
 975   par.forward=0;
 976   }*/
 977
 978   // Allocate memory for dynamic programming matrix
 979   hit.AllocateBacktraceMatrix(q.L+2,t.L+2); // ...with a separate dynamic programming matrix (memory!!)
 980   if (par.forward>=1 || Nstochali)
 981     hit.AllocateForwardMatrix(q.L+2,t.L+2);
 982   if (par.forward==2)
 983     hit.AllocateBackwardMatrix(q.L+2,t.L+2);
 984
 985 #ifndef CLUSTALO_NOFILE
 986   // Read structure file for Forward() function?
 987   if (strucfile && par.wstruc>0)
 988     {
 989       float PMIN=1E-20;
 990       Pstruc = new(float*[q.L+2]);
 991       for (int i=0; i<q.L+2; i++) Pstruc[i] = new(float[t.L+2]);
 992       Sstruc = new(float*[q.L+2]);
 993       for (int i=0; i<q.L+2; i++) Sstruc[i] = new(float[t.L+2]);
 994       FILE* strucf=NULL;
 995       if (strcmp(strucfile,"stdin"))
 996         {
 997           strucf = fopen(strucfile, "r");
 998           if (!strucf) OpenFileError(strucfile);
 999         }
1000       else
1001         {
1002           strucf = stdin;
1003           if (v>=2) printf("Reading structure matrix from standard input ... (for UNIX use ^D for 'end-of-file')\n");
1004         }
1005       for (int i=1; i<=q.L; i++)
1006         {
1007           for (int j=1; j<=t.L; j++)
1008             {
1009               float f;
1010               if (fscanf(strucf,"%f",&f) <=0 )
1011               {
1012                 fprintf(stderr,"Error: too few numbers in file %s while reading line %i, column %i\n",strucfile,i,j);
1013                 exit(1);
1014               }
1015               if (par.wstruc==1)
1016                 Pstruc[i][j]=fmax(f,PMIN);
1017               else
1018                 Pstruc[i][j]=fmax(pow(f,par.wstruc),PMIN);
1019 //            printf("%10.2E ",f);
1020               Sstruc[i][j] = par.wstruc * log2(f);
1021             }
1022 //        printf("\n");
1023         }
1024       fclose(strucf);
1025     } /* (strucfile && par.wstruc>0) */
1026 #endif
1027
1028   /* Do (self-)comparison, store results if score>SMIN,
1029      and try next best alignment */
1030   /* FIXME very ambigous and possibly faulty if-else */
1031   if (v>=2)
1032   {
1033     if (par.forward==2) {
1034         printf("Using maximum accuracy (MAC) alignment algorithm ...\n");
1035     }
1036     else if (par.forward==0) {
1037         printf("Using Viterbi algorithm ...\n");
1038     }
1039     else if (par.forward==1) {
1040         printf("Using stochastic sampling algorithm ...\n");
1041     }
1042     else {
1043         printf("\nWhat alignment algorithm are we using??\n");
1044     }
1045   }
1046
1047
1048   for (hit.irep=1; hit.irep<=imax(par.hitrank,par.altali); hit.irep++)
1049       {
1050           if (par.forward==0)
1051               {
1052                   // generate Viterbi alignment
1053                   hit.Viterbi(q,t,Sstruc);
1054                   hit.Backtrace(q,t);
1055               }
1056           else if (par.forward==1)
1057               {
1058                   // generate a single stochastically sampled alignment
1059                   hit.Forward(q,t,Pstruc);
1060                   srand( time(NULL) );   // initialize random generator
1061                   hit.StochasticBacktrace(q,t);
1062                   hitlist.Push(hit);      /* insert hit at beginning of list
1063                                              (last repeats first!) */
1064                   (hit.irep)++;
1065                   break;
1066               }
1067           else if (par.forward==2)
1068               {
1069                   // generate forward alignment
1070                   if (OK != hit.Forward(q,t,Pstruc)){
1071                       fprintf(stderr, "%s:%s:%d: cannot complete hit.Forward\n",
1072                               __FUNCTION__, __FILE__, __LINE__);
1073                       iRetVal = RETURN_FROM_MAC; /* spot double overflow in Forward(). FS, r241 -> r243 */
1074                       goto this_is_the_end;
1075                   }
1076                   if (OK != hit.Backward(q,t)){
1077                       fprintf(stderr, "%s:%s:%d: cannot complete hit.backward\n",
1078                               __FUNCTION__, __FILE__, __LINE__);
1079                       iRetVal = RETURN_FROM_MAC; /* spot double overflow in hit.Backward(). FS, r241 -> r243 */
1080                       goto this_is_the_end;
1081                   }
1082                   hit.MACAlignment(q,t);
1083                   if ((isnan(hit.score)) || (isnan(hit.Pforward))){
1084                       printf("nan after MAC\n");
1085                   }
1086                   hit.BacktraceMAC(q,t);
1087                   if ((isnan(hit.score)) || (isnan(hit.Pforward))){
1088                       printf("nan after backtrace\n");
1089                   }
1090               } /* use MAC algorithm */
1091           //       printf ("%-12.12s  %-12.12s   irep=%-2i  score=%6.2f hit.Pvalt=%.2g\n",hit.name,hit.fam,hit.irep,hit.score,hit.Pvalt);
1092           *dScore_p = hit.score;
1093
1094           if (hit.irep<=par.hitrank || hit.score>SMIN || (hit.Pvalt<pself && hit.score>0 ))
1095               hitlist.Push(hit);      // insert hit at beginning of list (last repeats first!) and do next alignment
1096           else
1097               {
1098                   if (hit.irep==1) hitlist.Push(hit); // first hit will be inserted into hitlist anyway, even if not significant
1099                   hit = hitlist.ReadLast(); // last alignment was not significant => read last (significant) hit from list
1100                   break;
1101                   /* FIXME: memory leak during normal alignment: both push'es above are not freed:
1102                    * valgrind --leak-check=full --show-reachable=yes ./src/clustalo -i ~/void/protein.fa -o /dev/null -v
1103                    ==9506== 1,456 bytes in 1 blocks are still reachable in loss record 4 of 5
1104                    ==9506==    at 0x4C2726C: operator new(unsigned long) (vg_replace_malloc.c:230)
1105                    ==9506==    by 0x4618C5: List<Hit>::Push(Hit) (list-C.h:134)
1106                    ==9506==    by 0x45E87A: hhalign (hhalign.cpp:914)
1107                    ==9506==    by 0x413F8C: HhalignWrapper (hhalign_wrapper.c:405)
1108                    ==9506==    by 0x41A7B5: MyMain (mymain.c:676)
1109                    ==9506==    by 0x4031A6: main (main.cpp:37)
1110                    ==9506==
1111                    ==9506==
1112                    ==9506== 4,442 (4,368 direct, 74 indirect) bytes in 3 blocks are definitely lost in loss record 5 of 5
1113                    ==9506==    at 0x4C2726C: operator new(unsigned long) (vg_replace_malloc.c:230)
1114                    ==9506==    by 0x4618C5: List<Hit>::Push(Hit) (list-C.h:134)
1115                    ==9506==    by 0x45E7C0: hhalign (hhalign.cpp:911)
1116                    ==9506==    by 0x413F8C: HhalignWrapper (hhalign_wrapper.c:405)
1117                    ==9506==    by 0x41A7B5: MyMain (mymain.c:676)
1118                    ==9506==    by 0x4031A6: main (main.cpp:37)
1119                   */
1120               }
1121       } /* 1 <= hit.irep <= imax(par.hitrank,par.altali) */
1122   Nali = hit.irep;
1123
1124
1125 #ifndef CLUSTALO_NOFILE
1126   // Write posterior probability matrix as TCoffee library file
1127   if (tcfile)
1128       {
1129           if (v>=2) printf("Writing TCoffee library file to %s\n",tcfile);
1130           int i,j;
1131           FILE* tcf=NULL;
1132           if (strcmp(tcfile,"stdout")) tcf = fopen(tcfile, "w"); else tcf = stdout;
1133           if (!tcf) OpenFileError(tcfile);
1134           fprintf(tcf,"! TC_LIB_FORMAT_01\n");
1135           fprintf(tcf,"%i\n",2); // two sequences in library file
1136           fprintf(tcf,"%s %i %s\n",q.name,q.L,q.seq[q.nfirst]+1);
1137           fprintf(tcf,"%s %i %s\n",hit.name,hit.L,hit.seq[hit.nfirst]+1);
1138           fprintf(tcf,"#1 2\n");
1139           for (i=1; i<=q.L; i++)  // print all pairs (i,j) with probability above PROBTCMIN
1140               for (j=1; j<=t.L; j++)
1141                   if (hit.B_MM[i][j]>probmin_tc)
1142                       fprintf(tcf,"%5i %5i %5i\n",i,j,iround(100.0*hit.B_MM[i][j]));
1143           for (int step=hit.nsteps; step>=1; step--)  // print all pairs on MAC alignment which were not yet printed
1144               {
1145                   i=hit.i[step]; j=hit.j[step];
1146                   //      printf("%5i %5i %5i  %i\n",i,j,iround(100.0*hit.B_MM[i][j]),hit.states[step]);
1147                   if (hit.states[step]>=MM && hit.B_MM[i][j]<=probmin_tc)
1148                       fprintf(tcf,"%5i %5i %5i\n",i,j,iround(100.0*hit.B_MM[i][j]));
1149               }
1150
1151
1152           fprintf(tcf,"! SEQ_1_TO_N\n");
1153           fclose(tcf);
1154           //       for (i=1; i<=q.L; i++)
1155           //            {
1156           //              double sum=0.0;
1157           //              for (j=1; j<=t.L; j++) sum+=hit.B_MM[i][j];
1158           //      printf("i=%-3i sum=%7.4f\n",i,sum);
1159           //            }
1160           //        printf("\n");
1161       } /* if (tcfile) */
1162
1163   // Write last alignment into alitabfile
1164   if (alitabfile)
1165       {
1166           FILE* alitabf=NULL;
1167           if (strcmp(alitabfile,"stdout")) alitabf = fopen(alitabfile, "w"); else alitabf = stdout;
1168           if (!alitabf) OpenFileError(alitabfile);
1169           if (par.forward==2)
1170               {
1171                   fprintf(alitabf,"    i     j  score     SS  probab\n");
1172                   for (int step=hit.nsteps; step>=1; step--)
1173                       if (hit.states[step]>=MM)
1174                           fprintf(alitabf,"%5i %5i %6.2f %6.2f %7.4f\n",hit.i[step],hit.j[step],hit.S[step],hit.S_ss[step],hit.P_posterior[step]);
1175               }
1176           else
1177               {
1178                   fprintf(alitabf,"    i     j  score     SS\n");
1179                   for (int step=hit.nsteps; step>=1; step--)
1180                       if (hit.states[step]>=MM)
1181                           fprintf(alitabf,"%5i %5i %6.2f %6.2f\n",hit.i[step],hit.j[step],hit.S[step],hit.S_ss[step]);
1182               }
1183           fclose(alitabf);
1184       } /* if (alitabfile) */
1185 #endif
1186
1187   // Do Stochastic backtracing?
1188   if (par.forward==1)
1189       for (int i=1; i<Nstochali; i++)
1190           {
1191               hit.StochasticBacktrace(q,t);
1192               hitlist.Push(hit);      //insert hit at beginning of list (last repeats first!)
1193               (hit.irep)++;
1194           }
1195   else // Set P-value, E-value and probability
1196       {
1197           if (q.mu) hitlist.GetPvalsFromCalibration(q);
1198           else if (t.mu) hitlist.GetPvalsFromCalibration(t);
1199       }
1200
1201   // Print FASTA or A2M alignments?
1202 #ifndef CLUSTALO_NOFILE
1203   if (*par.pairwisealisfile)
1204 #endif
1205       {
1206           if (v>=2) {
1207               cout<<"Printing alignments in " <<
1208                   (par.outformat==1? "FASTA" : par.outformat==2?"A2M" :"A3M") <<
1209                   " format to "<<par.pairwisealisfile<<"\n";
1210           }
1211           int iPrAliRtn = hitlist.PrintAlignments(
1212 #ifdef CLUSTALO
1213                                                   ppcFirstProf, ppcSecndProf, zcAux, zcError,
1214 #endif
1215                                                   q, par.pairwisealisfile, par.outformat);
1216           if ( (OK != iPrAliRtn) ){
1217               sprintf(zcAux, "%s:%s:%d: Could not print alignments\n",
1218                       __FUNCTION__, __FILE__, __LINE__);
1219               strcat(zcError, zcAux);
1220               iRetVal = RETURN_FROM_PRINT_ALI; /* this is where mis-alignment was originally spotted,
1221                                     hope to trap it now earlier. FS, r241 -> r243 */
1222               goto this_is_the_end;
1223           }
1224       } /* if (*par.pairwisealisfile) */
1225
1226 #ifndef CLUSTALO_NOFILE
1227   // Print hit list and  alignments
1228   if (*par.outfile)
1229       {
1230           hitlist.PrintHitList(q, par.outfile);
1231           hitlist.PrintAlignments(
1232 #ifdef CLUSTALO
1233                                   ppcFirstProf, ppcSecndProf, zcAux, zcError,
1234 #endif
1235                                   q, par.outfile);
1236           if (v>=2) {
1237               WriteToScreen(par.outfile,1000); //write only hit list to screen
1238           }
1239       }
1240 #endif
1241
1242
1243   ///////////////////////////////////////////////////////////////////////
1244 #ifndef CLUSTALO_NOFILE
1245   // Show results for hit with rank par.hitrank
1246   if (par.hitrank==0) hit=hitlist.Read(1); else hit=hitlist.Read(par.hitrank);
1247
1248   // Generate output alignment or HMM file?
1249   if (*par.alnfile || *par.psifile || *par.hhmfile)
1250       {
1251           if (par.append==0)
1252               {
1253                   if (v>=2 && *par.alnfile) printf("Merging template to query alignment and writing resulting alignment in A3M format to %s...\n",par.alnfile);
1254                   if (v>=2 && *par.psifile) printf("Merging template to query alignment and writing resulting alignment in PSI format to %s...\n",par.psifile);
1255               }
1256           else
1257               {
1258                   if (v>=2 && *par.alnfile) printf("Merging template to query alignment and appending template alignment in A3M format to %s...\n",par.alnfile);
1259                   if (v>=2 && *par.psifile) printf("Merging template to query alignment and appending template alignment in PSI format to %s...\n",par.psifile);
1260               }
1261
1262           // Read query alignment into Qali
1263           Alignment Qali;  // output A3M generated by merging A3M alignments for significant hits to the query alignment
1264           char qa3mfile[NAMELEN];
1265           RemoveExtension(qa3mfile,par.infile); // directory??
1266           strcat(qa3mfile,".a3m");
1267           FILE* qa3mf=fopen(qa3mfile,"r");
1268           if (!qa3mf) OpenFileError(qa3mfile);
1269           Qali.Read(qa3mf,qa3mfile);
1270           fclose(qa3mf);
1271
1272           // Align query with template in master-slave mode
1273           Qali.MergeMasterSlave(hit,par.tfile); /*FS, NOFILE2 (commented out) */
1274
1275           // Write output A3M alignment?
1276           if (*par.alnfile) {
1277               Qali.WriteToFile(par.alnfile,"a3m");
1278           }
1279
1280           if (*par.psifile)
1281               {
1282                   /* Convert ASCII to int (0-20),throw out all insert states,
1283                      record their number in I[k][i]  */
1284                   Qali.Compress("merged A3M file");
1285
1286                   // Write output PSI-BLAST-formatted alignment?
1287                   Qali.WriteToFile(par.psifile,"psi");
1288               }
1289       } /* if (*par.alnfile || *par.psifile || *par.hhmfile) */
1290 #endif  /* NOFILE2 */
1291   //////////////////////////////////////////////////////////////////////////
1292
1293
1294
1295
1296   //   double log2Pvalue;
1297   //   if (par.ssm && (par.ssm1 || par.ssm2))
1298   //     {
1299   //       log2Pvalue=hit.logPval/0.693147181+0.45*(4.0*par.ssw/0.15-hit.score_ss);
1300   //       if (v>=2)
1301   //    printf("Aligned %s with %s:\nApproximate P-value INCLUDING SS SCORE = %7.2g\n",q.name,t.name,pow(2.0,log2Pvalue));
1302   //     } else {
1303   //       if (v>=2)
1304   //    printf("Aligned %s with %s:\nApproximate P-value (without SS score) = %7.2g\n",q.name,t.name,hit.Pval);
1305   //    }
1306
1307   if (v>=2)
1308       {
1309           if (par.hitrank==0)
1310               printf("Aligned %s with %s: Score = %-7.2f  P-value = %-7.2g\n",
1311                      q.name,t.name,hit.score,hit.Pval);
1312           else
1313               printf("Aligned %s with %s (rank %i): Score = %-7.2f  P-value = %-7.2g\n",
1314                      q.name,t.name,par.hitrank,hit.score,hit.Pval);
1315       }
1316
1317   rHHscores->hhScore     = hit.score;
1318   rHHscores->forwardProb = hit.Pforward;
1319   rHHscores->sumPP       = hit.sum_of_probs;
1320
1321   /* next few lines commented out as they caused segfaults in RNA mode */
1322   /*  rHHscores->PP = (double *)malloc((hit.L+GOOD_MEASURE) * sizeof(double));*/
1323   rHHscores->L  = hit.L;
1324 /*  for (int i = 0; i < hit.L; i++){
1325       cout << "rHHscores->PP[i]" << rHHscores->PP[i] << endl;
1326       cout << "hit.P_posterior[i+1]" << hit.P_posterior[i+1] << endl;
1327       rHHscores->PP[i] = hit.P_posterior[i+1];
1328   } */
1329
1330
1331  this_is_the_end:
1332
1333   // Delete memory for dynamic programming matrix
1334   hit.DeleteBacktraceMatrix(q.L+2);
1335   if (par.forward>=1 || Nstochali)
1336       hit.DeleteForwardMatrix(q.L+2);
1337   if (par.forward==2)
1338       hit.DeleteBackwardMatrix(q.L+2);
1339   /*   if (Pstruc) {
1340        for (int i=0; i<q.L+2; i++) delete[](Pstruc[i]); delete[](Pstruc);}
1341   */
1342
1343   // Delete content of hits in hitlist
1344   hitlist.Reset();
1345   while (!hitlist.End())
1346       hitlist.ReadNext().Delete(); // Delete content of hit object
1347
1348   if (strucfile && par.wstruc>0)
1349       {
1350           for (int i=0; i<q.L+2; i++){
1351               delete[] Pstruc[i]; Pstruc[i] = NULL;
1352           }
1353           delete[] Pstruc; Pstruc = NULL;
1354           for (int i=0; i<q.L+2; i++){
1355               delete[] Sstruc[i]; Sstruc[i] = NULL;
1356           }
1357           delete[] Sstruc; Sstruc = NULL;
1358           delete[] strucfile; strucfile = NULL;
1359       }
1360
1361   if (pngfile){
1362       delete[] pngfile; pngfile = NULL;
1363   }
1364   if (alitabfile){
1365       delete[] alitabfile; alitabfile = NULL;
1366   }
1367   if (tcfile){
1368       delete[] tcfile; tcfile = NULL;
1369   }
1370   if (par.exclstr){
1371       delete[] par.exclstr; par.exclstr = NULL;
1372   }
1373
1374 #ifndef CLUSTALO_NOFILE
1375   // Print 'Done!'
1376   FILE* outf=NULL;
1377   if (!strcmp(par.outfile,"stdout")) printf("Done!\n");
1378   else
1379       {
1380           if (*par.outfile)
1381               {
1382                   outf=fopen(par.outfile,"a"); //open for append
1383                   fprintf(outf,"Done!\n");
1384                   fclose(outf);
1385               }
1386           if (v>=2) printf("Done\n");
1387       }
1388 #endif
1389
1390   //qali.ClobberGlobal();
1391   hit.ClobberGlobal();
1392   if (iFirstCnt > 0){
1393       q.ClobberGlobal();
1394   }
1395   if (iSecndCnt > 0){
1396       t.ClobberGlobal();
1397   }
1398   hitlist.ClobberGlobal();
1399
1400   return iRetVal;
1401
1402 } /* this is the end of hhalign() //end main */
1403
1404 //////////////////////////////////////////////////////////////////////////////
1405 // END OF MAIN
1406 //////////////////////////////////////////////////////////////////////////////
1407
1408 /*
1409  * EOF hhalign.C
1410  */
1411