@@ -77,15 +77,15 @@ Alignment::Alignment(int maxseq, int maxres)

     //printf(">>>>>>>>%s:%s:%d: maxseq=%d, maxres=%d\n", __FUNCTION__, __FILE__, __LINE__, maxseq, maxres); /* (FS) */

   longname = new(char[DESCLEN]);

-  sname = new(char*[maxseq+2]); /* MR1 */

-  seq = new(char*[maxseq+2]); /* MR1 */

-  l = new(int[maxres]);

-  X = new(char*[maxseq+2]);  /* MR1 */

-  I = new(short unsigned int*[maxseq+2]); /* MR1 */

-  keep = new(char[maxseq+2]); /* MR1 */

-  display = new(char[maxseq+2]); /* MR1 */

-  wg = new(float[maxseq+2]); /* MR1 */

-  nseqs = new(int[maxres+2]); /* MR1 */

+  sname = new char*[maxseq+2]; /* MR1 */

+  seq = new char*[maxseq+2]; /* MR1 */

+  l = new int[maxres];

+  X = new char*[maxseq+2];  /* MR1 */

+  I = new short unsigned int*[maxseq+2]; /* MR1 */

+  keep = new char[maxseq+2]; /* MR1 */

+  display = new char[maxseq+2]; /* MR1 */

+  wg = new float[maxseq+2]; /* MR1 */

+  nseqs = new int[maxres+2]; /* MR1 */

   N_in=L=0;

   nres=NULL; // number of residues per sequence k

   first=NULL; // first residue in sequence k

@@ -138,7 +138,7 @@ Alignment::Read(FILE* inf, char infile[], char* firstline)

   int k; // Index of sequence being read currently (first=0)

   char line[LINELEN]=""; // input line

   //char cur_seq[MAXCOL]; // Sequence currently read in

-  char *cur_seq=new(char[par.maxColCnt]);

+  char *cur_seq=new char[par.maxColCnt];

   char* cur_name; // Sequence currently read in

   int linenr=0; // current line number in input file

   char skip_sequence=0;

@@ -180,11 +180,11 @@ Alignment::Read(FILE* inf, char infile[], char* firstline)

                               }

                           // Create space for residues and paste new sequence in

-                          seq[k]=new(char[strlen(cur_seq)+2]);

+                          seq[k]=new char[strlen(cur_seq)+2];

                           if (!seq[k]) MemoryError("array for input sequences");

-                          X[k]=new(char[strlen(cur_seq)+2]);

+                          X[k]=new char[strlen(cur_seq)+2];

                           if (!X[k]) MemoryError("array for input sequences");

-                          I[k]=new(short unsigned int[strlen(cur_seq)+2]);

+                          I[k]=new short unsigned int[strlen(cur_seq)+2];

                           if (!I[k]) MemoryError("array for input sequences");

                           strcpy(seq[k],cur_seq);

                       }

@@ -233,7 +233,7 @@ Alignment::Read(FILE* inf, char infile[], char* firstline)

                   // store sequence name

                   if (v>=4) printf("Reading seq %-16.16s k=%3i n_displ=%3i display[k]=%i keep[k]=%i\n",cur_name,k,n_display,display[k],keep[k]);

-                  sname[k] = new(char[strlen(cur_name)+1]);

+                  sname[k] = new char[strlen(cur_name)+1];

                   if (!sname[k]) {MemoryError("array for sequence names");}

                   strcpy(sname[k],cur_name);

               } // end if(line contains sequence name)

@@ -348,11 +348,11 @@ Alignment::Read(FILE* inf, char infile[], char* firstline)

   if (k>=0) //if at least one sequence was read in

       {

-          seq[k]=new(char[strlen(cur_seq)+2]);

+          seq[k]=new char[strlen(cur_seq)+2];

           if (!seq[k]) MemoryError("array for input sequences");

-          X[k]=new(char[strlen(cur_seq)+2]);

+          X[k]=new char[strlen(cur_seq)+2];

           if (!X[k]) MemoryError("array for input sequences");

-          I[k]=new(short unsigned int[strlen(cur_seq)+2]);

+          I[k]=new short unsigned int[strlen(cur_seq)+2];

           if (!I[k]) MemoryError("array for input sequences");

           strcpy(seq[k],cur_seq);

       }

@@ -420,7 +420,7 @@ Alignment::Compress(const char infile[])

     /*static short unsigned int h[MAXSEQ];*/

     /*short*/ unsigned int *h = NULL; /* short may lead to overflow for long alignments, FS, r235 -> r236 */

-    h = new(/*short*/ unsigned int[N_in+2]); /* short -> overflow, FS, r235 -> r236 */

+    h = new /*short*/ unsigned int[N_in+2]; /* short -> overflow, FS, r235 -> r236 */

     float *percent_gaps = NULL; /* FS, 2010-Nov */

     char *match_state = NULL;  /* FS, 2010-Nov */

@@ -552,13 +552,13 @@ Alignment::Compress(const char infile[])

                had to move declaration of float *percent_gaps out of switch()

             */

             //float percent_gaps[MAXCOL]; //percentage of gaps in column k (with weighted sequences)

-            percent_gaps = new(float[par.maxColCnt]);

+            percent_gaps = new float[par.maxColCnt];

             //determine number of columns L in alignment

             L=strlen(seq[kfirst])-1;

             // Conversion to integer representation, checking for unequal lengths and initialization

-            if (nres==NULL) nres=new(int[N_in]);

+            if (nres==NULL) nres=new int[N_in];

             for (k=0; k<N_in; k++)

                 {

                     if (!keep[k]) continue;

@@ -780,7 +780,7 @@ Alignment::Compress(const char infile[])

                had to move declaration of float *percent_gaps out of switch()

             */

             //char match_state[MAXCOL]; //1: column assigned to match state 0: insert state

-            match_state = new(char[par.maxColCnt]);

+            match_state = new char[par.maxColCnt];

             // Determine number of columns L in alignment

             L=strlen(seq[0]+1);

@@ -942,7 +942,7 @@ Alignment::FilterForDisplay(int max_seqid, int coverage, int qid, float qsc, int

     if (par.mark) return n_display;

-    char *dummy = new(char[N_in+1]);

+    char *dummy = new char[N_in+1];

     int vtmp=v, seqid;

     v=0;

     n_display=0;

@@ -1005,12 +1005,12 @@ Alignment::Filter2(char keep[], int coverage, int qid, float qsc, int seqid1, in

 {

     // In the beginnning, keep[k] is 1 for all regular amino acid sequences and 0 for all others (ss_conf, ss_pred,...)

     // In the end, keep[k] will be 1 for all regular representative sequences kept in the alignment, 0 for all others

-    char* in=new(char[N_in+1]); // in[k]=1: seq k has been accepted; in[k]=0: seq k has not yet been accepted at current seqid

-    char* inkk=new(char[N_in+1]); // inkk[k]=1 iff in[ksort[k]]=1 else 0;

-    int* Nmax=new(int[L+2]); // position-dependent maximum-sequence-identity threshold for filtering /* MR1, used to be called idmax*/

-    int* idmaxwin=new(int[L+2]); // minimum value of idmax[i-WFIL,i+WFIL]

-    int* seqid_prev=new(int[N_in+1]); // maximum-sequence-identity threshold used in previous round of filtering (with lower seqid)

-    int* N=new(int[L+2]); // N[i] number of already accepted sequences at position i

+    char* in=new char[N_in+1]; // in[k]=1: seq k has been accepted; in[k]=0: seq k has not yet been accepted at current seqid

+    char* inkk=new char[N_in+1]; // inkk[k]=1 iff in[ksort[k]]=1 else 0;

+    int* Nmax=new int[L+2]; // position-dependent maximum-sequence-identity threshold for filtering /* MR1, used to be called idmax*/

+    int* idmaxwin=new int[L+2]; // minimum value of idmax[i-WFIL,i+WFIL]

+    int* seqid_prev=new int[N_in+1]; // maximum-sequence-identity threshold used in previous round of filtering (with lower seqid)

+    int* N=new int[L+2]; // N[i] number of already accepted sequences at position i

     const int WFIL=25; // see previous line

     int diffNmax=Ndiff;       // current  maximum difference of Nmax[i] and Ndiff /* MR1 */

@@ -1039,8 +1039,8 @@ Alignment::Filter2(char keep[], int coverage, int qid, float qsc, int seqid1, in

     // Determine first[k], last[k]?

     if (first==NULL)

         {

-            first=new(int[N_in]);// first non-gap position in sequence k

-            last =new(int[N_in]);// last  non-gap position in sequence k

+            first=new int[N_in];// first non-gap position in sequence k

+            last =new int[N_in];// last  non-gap position in sequence k

             for (k=0; k<N_in; k++) // do this for ALL sequences, not only those with in[k]==1 (since in[k] may be display[k])

                 {

                     for (i=1; i<=L; i++) if (X[k][i]<NAA) break;

@@ -1054,7 +1054,7 @@ Alignment::Filter2(char keep[], int coverage, int qid, float qsc, int seqid1, in

     if ( (nres==NULL)  || (sizeof(nres)<N_in*sizeof(int)) )

         {

     	 if (nres) delete[] nres;

-    		nres=new(int[N_in]);

+    		nres=new int[N_in];

             for (k=0; k<N_in; k++) // do this for ALL sequences, not only those with in[k]==1 (since in[k] may be display[k])

                 {

                     int nr=0;

@@ -1068,7 +1068,7 @@ Alignment::Filter2(char keep[], int coverage, int qid, float qsc, int seqid1, in

     // Sort sequences according to length; afterwards, nres[ksort[kk]] is sorted by size

     if (ksort==NULL)

         {

-            ksort=new(int[N_in]); // never reuse alignment object for new alignment with more sequences

+            ksort=new int[N_in]; // never reuse alignment object for new alignment with more sequences

             for (k=0; k<N_in; k++) ksort[k]=k;

             QSortInt(nres,ksort,kfirst+1,N_in-1,-1); //Sort sequences after kfirst (query) in descending order

         }

@@ -1311,7 +1311,7 @@ Alignment::HomologyFilter(int coverage_core, float qsc_core, float coresc)

     const int Ndiff_core=0;

     int n;

     HMM qcore;

-    char* coreseq=new(char[N_in]); // coreseq[k]=1 if sequence belongs to core of alignment (i.e. it is very similar to query)

+    char* coreseq=new char[N_in]; // coreseq[k]=1 if sequence belongs to core of alignment (i.e. it is very similar to query)

     for (int k=0; k<N_in; k++) coreseq[k]=keep[k]; // Copy keep[] into coreseq[]

     // Remove sequences with seq. identity larger than seqid percent (remove the shorter of two)

@@ -1362,7 +1362,7 @@ Alignment::FilterWithCoreHMM(char in[], float coresc, HMM& qcore)

     int i; // column in query alignment

     int a; // amino acid (0..19)

     int n=1; // number of sequences that passed filter

-    float** logodds=new(float*[L+1]); // log-odds ratios for HMM qcore

+    float** logodds=new float*[L+1]; // log-odds ratios for HMM qcore

     char gap; // 1: previous state in seq k was a gap 0: previous state in seq k was an amino acid

     float score; // score of sequence k aligned with qcore

@@ -1371,8 +1371,8 @@ Alignment::FilterWithCoreHMM(char in[], float coresc, HMM& qcore)

     // Determine first[k], last[k]?

     if (first==NULL)

         {

-            first=new(int[N_in]);// first non-gap position in sequence k

-            last =new(int[N_in]);// last non-gap position in sequence k

+            first=new int[N_in];// first non-gap position in sequence k

+            last =new int[N_in];// last non-gap position in sequence k

             for (k=0; k<N_in; k++) // do this for ALL sequences, not only those with in[k]==1 (since in[k] may be display[k])

                 {

                     for (i=1; i<=L; i++) if (X[k][i]<NAA) break;

@@ -1385,7 +1385,7 @@ Alignment::FilterWithCoreHMM(char in[], float coresc, HMM& qcore)

     // Determine number of residues nres[k]?

     if (nres==NULL)

         {

-            nres=new(int[N_in]);

+            nres=new int[N_in];

             for (k=0; k<N_in; k++) // do this for ALL sequences, not only those with in[k]==1 (since in[k] may be display[k])

                 {

                     int nr=0;

@@ -1667,18 +1667,18 @@ Alignment::FrequenciesAndTransitions(HMM& q, char* in)

                     q.sname[q.ncons]=new(char[10]);

                     if (!q.sname[q.ncons]) {MemoryError("array of names for displayed sequences");}

                     strcpy(q.sname[q.ncons],"Consensus");

-                    q.seq[q.ncons]=new(char[L+2]);

+                    q.seq[q.ncons]=new char[L+2];

                     if (!q.seq[q.ncons]) {MemoryError("array of names for displayed sequences");}

                 }

             if (par.cons)

                 {

                     // Reserve space for consensus sequence as first sequence in alignment

                     q.nfirst=n++; kfirst=-1;

-                    q.sname[q.nfirst]=new(char[strlen(name)+11]);

+                    q.sname[q.nfirst]=new char[strlen(name)+11];

                     if (!q.sname[q.nfirst]) {MemoryError("array of names for displayed sequences");}

                     strcpy(q.sname[q.nfirst],name);

                     strcat(q.sname[q.nfirst],"_consensus");

-                    q.seq[q.nfirst]=new(char[L+2]);

+                    q.seq[q.nfirst]=new char[L+2];

                     if (!q.seq[q.nfirst]) {MemoryError("array of names for displayed sequences");}

                 }

             // Calculate consensus amino acids using similarity matrix

@@ -1736,10 +1736,10 @@ Alignment::FrequenciesAndTransitions(HMM& q, char* in)

                     else if (k==kfirst) nn=q.nfirst=n++;

                     else nn=n++;

                     // strcut(sname[k]," "); // delete rest of name line beginning with two spaces " " // Why this?? Problem for pdb seqs without chain

-                    q.sname[nn]=new(char[strlen(sname[k])+1]);

+                    q.sname[nn]=new char[strlen(sname[k])+1];

                     if (!q.sname[nn]) {MemoryError("array of names for displayed sequences");}

                     strcpy(q.sname[nn],sname[k]);

-                    q.seq[nn]=new(char[strlen(seq[k])+1]);

+                    q.seq[nn]=new char[strlen(seq[k])+1];

                     if (!q.seq[nn]) {MemoryError("array of names for displayed sequences");}

                     strcpy(q.seq[nn],seq[k]);

                 }

@@ -1826,14 +1826,14 @@ Alignment::Amino_acid_frequencies_and_transitions_from_M_state(HMM& q, char* in)

   //float wi[MAXSEQ]; // weight of sequence k in column i, calculated from subalignment i

   float *wi=NULL; // weight of sequence k in column i, calculated from subalignment i

   //float Neff[MAXRES]; // diversity of subalignment i

-  float *Neff = new(float[par.maxResLen]); // diversity of subalignment i

+  float *Neff = new float[par.maxResLen]; // diversity of subalignment i

   int nseqi=0; // number of sequences in subalignment i

   int ncol=0; // number of columns j that contribute to Neff[i]

   char change; // has the set of sequences in subalignment changed? 0:no 1:yes

   float fj[NAA+3]; // to calculate entropy

   float sum;

-  wi = new(float[N_in+2]);

+  wi = new float[N_in+2];

   // Global weights?

   if (par.wg==1)

@@ -1842,7 +1842,7 @@ Alignment::Amino_acid_frequencies_and_transitions_from_M_state(HMM& q, char* in)

   // Initialization

   q.Neff_HMM=0.0f;

   Neff[0]=0.0; // if the first column has no residues (i.e. change==0), Neff[i]=Neff[i-1]=Neff[0]

-  n = new(int*[L+2]);

+  n = new int*[L+2];

   for (j=1; j<=L; j++) n[j]=new(int[NAA+3]);

   for (j=1; j<=L; j++)

     for (a=0; a<NAA+3; a++) n[j][a]=0;

@@ -2052,7 +2052,7 @@ Alignment::Transitions_from_I_state(HMM& q, char* in)

     //float wi[MAXSEQ]; // weight of sequence k in column i, calculated from subalignment i

     float *wi = NULL; // weight of sequence k in column i, calculated from subalignment i

     //float Neff[MAXRES]; // diversity of subalignment i

-    float *Neff = new(float[par.maxResLen]); // diversity of subalignment i

+    float *Neff = new float[par.maxResLen]; // diversity of subalignment i

     int nseqi; // number of sequences in subalignment i

     int ncol; // number of columns j that contribute to Neff[i]

     float fj[NAA+3]; // to calculate entropy

@@ -2060,7 +2060,7 @@ Alignment::Transitions_from_I_state(HMM& q, char* in)

     float Nlim=0.0; // only for global weights

     float scale=0.0; // only for global weights

-    wi = new(float[N_in+2]);

+    wi = new float[N_in+2];

     // Global weights?

     if (par.wg==1)

@@ -2071,7 +2071,7 @@ Alignment::Transitions_from_I_state(HMM& q, char* in)

         }

     // Initialization

-    n = new(int*[L+2]);

+    n = new int*[L+2];

     for (j=1; j<=L; j++) n[j]=new(int[NAA+3]);

     //////////////////////////////////////////////////////////////////////////////////////////////

@@ -2254,7 +2254,7 @@ Alignment::Transitions_from_D_state(HMM& q, char* in)

     //float wi[MAXSEQ]; // weight of sequence k in column i, calculated from subalignment i

     float *wi=NULL; // weight of sequence k in column i, calculated from subalignment i

     //float Neff[MAXRES]; // diversity of subalignment i 

-    float *Neff = new(float[par.maxResLen]); // diversity of subalignment i 

+    float *Neff = new float[par.maxResLen]; // diversity of subalignment i

     int nseqi=0; // number of sequences in subalignment i (for DEBUGGING)

     int ncol=0; // number of columns j that contribute to Neff[i]

     char change; // has the set of sequences in subalignment changed? 0:no 1:yes

@@ -2263,7 +2263,7 @@ Alignment::Transitions_from_D_state(HMM& q, char* in)

     float Nlim=0.0; // only for global weights

     float scale=0.0; // only for global weights

-    wi = new(float[N_in+2]); /* FIXME: FS */

+    wi = new float[N_in+2]; /* FIXME: FS */

     // Global weights?

     if (par.wg==1)

         {

@@ -2273,7 +2273,7 @@ Alignment::Transitions_from_D_state(HMM& q, char* in)

         }

     // Initialization

-    n = new(int*[L+2]);

+    n = new int*[L+2];

     for (j=1; j<=L; j++) n[j]=new(int[NAA+3]);

     for (j=1; j<=L; j++)

         for (a=0; a<NAA+3; a++) n[j][a]=0;

@@ -2544,7 +2544,7 @@ Alignment::MergeMasterSlave(Hit& hit, char ta3mfile[])

  N_filtered = Tali.Filter(par.max_seqid,par.coverage,par.qid,par.qsc,par.Ndiff);

  // Record imatch[j]

- int* imatch=new(int[hit.j2+1]);

+ int* imatch=new int[hit.j2+1];

  int step = hit.nsteps;

  for (j=hit.j1; j<=hit.j2; j++)

  {

@@ -2676,7 +2676,7 @@ Alignment::MergeMasterSlave(Hit& hit, char ta3mfile[])

  iprev=i; lprev=l;

  if (h>=maxcol-1000) // too few columns? Reserve double space

  {

- char* new_seq=new(char[2*maxcol]);

+ char* new_seq=new char[2*maxcol];

  strncpy(new_seq,cur_seq,maxcol); //////// check: maxcol-1 ????

  delete[](cur_seq); (cur_seq) = NULL;

  cur_seq=new_seq;

@@ -2689,14 +2689,14 @@ Alignment::MergeMasterSlave(Hit& hit, char ta3mfile[])

  cur_seq[h++]='\0';

  keep[N_in] = display[N_in] = KEEP_CONDITIONALLY;

- seq[N_in]=new(char[h]);

+ seq[N_in]=new char[h];

  if (!seq[N_in]) MemoryError("array for input sequences");

  strcpy(seq[N_in],cur_seq);

- X[N_in]=new(char[h]);

+ X[N_in]=new char[h];

  if (!X[N_in]) MemoryError("array for input sequences");

- I[N_in]=new(short unsigned int[h]);

+ I[N_in]=new short unsigned int[h];

  if (!I[N_in]) MemoryError("array for input sequences");

- sname[N_in]=new(char[strlen(Tali.sname[k])+1]);

+ sname[N_in]=new char[strlen(Tali.sname[k])+1];

  if (!sname[N_in]) MemoryError("array for input sequences");

  strcpy(sname[N_in],Tali.sname[k]);

  N_in++;

@@ -2727,7 +2727,7 @@ Alignment::AddSequence(char Xk[], int Ik[])

 {

  int i; // position in query and target

  if (L<=0) InternalError("L is not set in AddSequence()");

- X[N_in]=new(char[L+2]);

+ X[N_in]=new char[L+2];

  for (i=0; i<=L+1; i++) X[N_in][i]=Xk[i];

  if (Ik==NULL)

  for (i=0; i<=L+1; i++) I[N_in][i]=0;

@@ -2777,7 +2777,7 @@ Alignment::GetPositionSpecificWeights(float* w[])

  {

  // Initialization

- n = new(int*[L+2]);

+ n = new int*[L+2];

  for (j=1; j<=L; j++) n[j]=new(int[NAA+3]);

  for (j=1; j<=L; j++)

  for (a=0; a<NAA+3; a++) n[j][a]=0;

@@ -2894,9 +2894,9 @@ Alignment::Transfer(char **ppcProf, int iCnt){

     /* @<allocate memory for sequences etc@> */

     for (k = 0; k < iCnt; k++){

 #define GOOD_MEASURE 1000 /* Temporary -- can be removed once rest in place */

-        I[k] = new(short unsigned int[iLen+2+GOOD_MEASURE]);

-        X[k] = new(char[iLen+2+GOOD_MEASURE]);

-        seq[k] = new(char[iLen+2+GOOD_MEASURE]);

+        I[k] = new short unsigned int[iLen+2+GOOD_MEASURE];

+        X[k] = new char[iLen+2+GOOD_MEASURE];

+        seq[k] = new char[iLen+2+GOOD_MEASURE];

         seq[k][0] = ' ';

         seq[k][1] = '\0';

         if (NULL == ppcProf[k]){

new file mode 100644

@@ -0,0 +1,3029 @@

+/* -*- mode: c; tab-width: 4; c-basic-offset: 4; indent-tabs-mode: nil -*- */

+

+/*********************************************************************

+ * Clustal Omega - Multiple sequence alignment

+ *

+ * Copyright (C) 2010 University College Dublin

+ *

+ * Clustal-Omega is free software; you can redistribute it and/or

+ * modify it under the terms of the GNU General Public License as

+ * published by the Free Software Foundation; either version 2 of the

+ * License, or (at your option) any later version.

+ *

+ * This file is part of Clustal-Omega.

+ *

+ ********************************************************************/

+

+/*

+ *  RCS $Id: hhalignment-C.h 236 2011-04-14 11:30:04Z fabian $

+ */

+

+

+/*

+ * Changelog: Michael Remmert made changes to hhalign stand-alone code

+ * FS implemented some of the changes on 2010-10-28 -> MR1

+ *

+ * Note: MR seems to have changed all [aijk]++ to ++[aijk],

+ * FS did not do that on 2010-10-28

+ */

+

+// hhalignment.C

+

+//////////////////////////////////////////////////////////////////////////////

+//// Class Alignment

+//////////////////////////////////////////////////////////////////////////////

+

+// hhalignment.C

+

+#ifndef MAIN

+#define MAIN

+#include <iostream> // cin, cout, cerr

+#include <fstream> // ofstream, ifstream

+#include <stdio.h> // printf

+using std::cout;

+using std::cerr;

+using std::endl;

+using std::ios;

+using std::ifstream;

+using std::ofstream;

+#include <stdlib.h> // exit

+#include <string> // strcmp, strstr

+#include <math.h> // sqrt, pow

+#include <limits.h> // INT_MIN

+#include <float.h> // FLT_MIN

+#include <time.h> // clock

+#include <ctype.h> // islower, isdigit etc

+#include "util-C.h" // imax, fmax, iround, iceil, ifloor, strint, strscn, strcut, substr, uprstr, uprchr, Basename etc.

+#include "list.h" // list data structure

+#include "hash.h" // hash data structure

+#include "hhdecl-C.h"

+#include "hhutil-C.h" // imax, fmax, iround, iceil, ifloor, strint, strscn, strcut, substr, uprstr, uprchr, Basename etc.

+#include "hhhmm.h"

+#endif

+

+

+enum {KEEP_NOT = 0, KEEP_CONDITIONALLY, KEEP_ALWAYS};

+

+//////////////////////////////////////////////////////////////////////////////

+// Class Alignment

+//////////////////////////////////////////////////////////////////////////////

+

+

+//////////////////////////////////////////////////////////////////////////////

+// Object constructor

+//////////////////////////////////////////////////////////////////////////////

+Alignment::Alignment(int maxseq, int maxres)

+{

+

+    //printf(">>>>>>>>%s:%s:%d: maxseq=%d, maxres=%d\n", __FUNCTION__, __FILE__, __LINE__, maxseq, maxres); /* (FS) */

+  longname = new(char[DESCLEN]);

+  sname = new(char*[maxseq+2]); /* MR1 */

+  seq = new(char*[maxseq+2]); /* MR1 */

+  l = new(int[maxres]);

+  X = new(char*[maxseq+2]);  /* MR1 */

+  I = new(short unsigned int*[maxseq+2]); /* MR1 */

+  keep = new(char[maxseq+2]); /* MR1 */

+  display = new(char[maxseq+2]); /* MR1 */

+  wg = new(float[maxseq+2]); /* MR1 */

+  nseqs = new(int[maxres+2]); /* MR1 */

+  N_in=L=0;

+  nres=NULL; // number of residues per sequence k

+  first=NULL; // first residue in sequence k

+  last=NULL; // last residue in sequence k

+  ksort=NULL; // sequence indices sorted by descending nres[k]

+  name[0]='\0'; // no name defined yet

+  longname[0]='\0'; // no name defined yet

+  fam[0]='\0'; // no name defined yet

+  file[0]='\0'; // no name defined yet

+  readCommentLine = '0'; /* MR1 */

+}

+

+//////////////////////////////////////////////////////////////////////////////

+// Object destructor

+//////////////////////////////////////////////////////////////////////////////

+Alignment::~Alignment()

+{

+  delete[] longname; longname = NULL;

+  for(int k=0; k<N_in; k++)

+    {

+      delete[] sname[k]; sname[k] = NULL;

+      delete[] seq[k]; seq[k] = NULL;

+      delete[] X[k]; X[k] = NULL;

+      delete[] I[k]; I[k] = NULL;

+    }

+  delete[] sname; sname = NULL;

+  delete[] seq; seq = NULL;

+  delete[] X; X = NULL;

+  delete[] I; I = NULL;

+  delete[] l; l = NULL;

+  delete[] keep; keep = NULL;

+  delete[] display; display = NULL;

+  delete[] wg; wg = NULL;

+  delete[] nseqs; nseqs = NULL;

+  delete[] nres; nres = NULL;

+  delete[] first; first = NULL;

+  delete[] last; last = NULL;

+  delete[] ksort; ksort = NULL;

+}

+

+

+/**

+ * @brief Reads in an alignment from file into matrix seq[k][l] as ASCII

+ */

+void 

+Alignment::Read(FILE* inf, char infile[], char* firstline)

+{

+  int l; // Postion in alignment incl. gaps (first=1)

+  int h; // Position in input line (first=0)

+  int k; // Index of sequence being read currently (first=0)

+  char line[LINELEN]=""; // input line

+  //char cur_seq[MAXCOL]; // Sequence currently read in

+  char *cur_seq=new(char[par.maxColCnt]);

+  char* cur_name; // Sequence currently read in

+  int linenr=0; // current line number in input file

+  char skip_sequence=0;

+  RemoveExtension(file,infile); //copy rootname (w/o path) of infile into file variable of class object

+

+  kss_dssp=ksa_dssp=kss_pred=kss_conf=kfirst=-1;

+  n_display=0;

+  N_in=0;

+  N_filtered=0;

+  N_ss=0;

+  cur_seq[0]=' '; // overwrite '\0' character at beginning to be able to do strcpy(*,cur_seq)

+  l=1; k=-1;

+

+  // Does firstline already contain first line of file?

+  if (firstline!= NULL) strcpy(line,firstline);

+

+  /////////////////////////////////////////////////////////////////////////

+  // Read infile line by line

+  /* FIXME: not safe to use MAXSEQ, however, don't think we ever get here (FS) */

+  while(firstline || (fgetline(line,LINELEN,inf) && (k<MAXSEQ))) /* FIXME: FS introduced () around &&, precedence! MR1 */

+      {

+          linenr++;

+          firstline=NULL;

+          if (line[0]=='>') //line contains sequence name

+              {

+                  if (k>=MAXSEQ-1)

+                      {

+                          if (v>=1 && k>=MAXSEQ)

+                              cerr<<endl<<"WARNING: maximum number "<<MAXSEQ<<" of sequences exceded in file "<<infile<<"\n";

+                          break;

+                      }

+                  cur_name=line+1; //beginning of current sequence name

+                  if (k>=0) //if this is at least the second name line

+                      {

+                          if (strlen(cur_seq)==0)

+                              {

+                                  cerr<<endl<<"Error: sequence "<<sname[k]<<" contains no residues."<<endl;

+                                  throw 1;

+                              }

+

+                          // Create space for residues and paste new sequence in

+                          seq[k]=new(char[strlen(cur_seq)+2]);

+                          if (!seq[k]) MemoryError("array for input sequences");

+                          X[k]=new(char[strlen(cur_seq)+2]);

+                          if (!X[k]) MemoryError("array for input sequences");

+                          I[k]=new(short unsigned int[strlen(cur_seq)+2]);

+                          if (!I[k]) MemoryError("array for input sequences");

+                          strcpy(seq[k],cur_seq);

+                      }

+                  skip_sequence=0;

+

+                  k++;

+                  l=1; //position in current sequence (first=1)

+

+                  // display[k]= 0: do not show in Q-T alignments 1: show if not filtered out later 2: show in any case (do not filter out)

+                  // keep[k] = 0: do not include in profile 1: include if not filtered out later 2: include in any case (do not filter out)

+                  /* {KEEP_NOT=0, KEEP_CONDITIONALLY=1, KEEP_ALWAYS=2} */

+                  if (line[1]=='@') cur_name++; //skip @-character in name

+                  if (!strncmp(line,">ss_dssp",8)) {

+                      if (kss_dssp<0) {display[k]=2; n_display++; keep[k]=KEEP_NOT; kss_dssp=k; N_ss++;} else {skip_sequence=1; k--; continue;}

+                  }

+                  else if (!strncmp(line,">sa_dssp",8)) {

+                      if (ksa_dssp<0) {display[k]=KEEP_ALWAYS; n_display++; keep[k]=KEEP_NOT; ksa_dssp=k; N_ss++;} else {skip_sequence=1; k--; continue;}

+                  }

+                  else if (!strncmp(line,">ss_pred",8)) {

+                      if (kss_pred<0) {display[k]=KEEP_ALWAYS; n_display++; keep[k]=KEEP_NOT; kss_pred=k; N_ss++;} else {skip_sequence=1; k--; continue;}

+                  }

+                  else if (!strncmp(line,">ss_conf",8)) {

+                      if (kss_conf<0) {display[k]=KEEP_ALWAYS; n_display++; keep[k]=KEEP_NOT; kss_conf=k; N_ss++;} else {skip_sequence=1; k--; continue;}

+                  }

+                  else if (!strncmp(line,">ss_",4) || !strncmp(line,">sa_",4)) {

+                      display[k]=KEEP_ALWAYS; n_display++; keep[k]=KEEP_NOT; N_ss++;

+                  }

+                  else if (!strncmp(line,">aa_",4)) { // ignore sequences beginning with ">aa_"

+                      skip_sequence=1; k--; continue;

+                  }

+                  //store first real seq

+                  else if (kfirst<0)

+                      {

+                          char word[NAMELEN];

+                          strwrd(word,line); // Copies first word in ptr to str

+                          if (strstr(word,"_consensus"))

+                              {display[k]=2; keep[k]=0; n_display++; kfirst=k;} /* MR1 */

+                          else

+                              {display[k]=keep[k]=KEEP_ALWAYS; n_display++; kfirst=k;}

+                      }

+                  //store all sequences

+                  else if (par.mark==0) {display[k]=keep[k]=KEEP_CONDITIONALLY; n_display++;}

+                  //store sequences up to nseqdis

+                  else if (line[1]=='@'&& n_display-N_ss<par.nseqdis) {display[k]=keep[k]=KEEP_ALWAYS; n_display++;}

+                  else {display[k]=KEEP_NOT; keep[k]=KEEP_CONDITIONALLY;}

+

+                  // store sequence name

+                  if (v>=4) printf("Reading seq %-16.16s k=%3i n_displ=%3i display[k]=%i keep[k]=%i\n",cur_name,k,n_display,display[k],keep[k]);

+                  sname[k] = new(char[strlen(cur_name)+1]);

+                  if (!sname[k]) {MemoryError("array for sequence names");}

+                  strcpy(sname[k],cur_name);

+              } // end if(line contains sequence name)

+

+          else if (line[0]=='#') // Commentary line?

+              {

+                  // #PF01367.9 5_3_exonuc: 5'-3' exonuclease, C-terminal SAM fold; PDB 1taq, 1bgx (T:271-174), 1taq (271-174)

+                  if (name[0]) continue; // if already name defined: skip commentary line

+                  char *ptr1, *ptr2;

+                  ptr1=strscn_(line+1); // set ptr1 to first non-whitespace character after '#' -> AC number

+                  strncpy(longname,ptr1,DESCLEN-1); // copy whole commentary line after '# ' into longname

+                  longname[DESCLEN-1]='\0';

+                  strtr(longname,""," ");

+                  ptr2=strcut_(ptr1); // cut after AC number and set ptr2 to first non-whitespace character after AC number

+                  // strcpy(fam,ptr1); // copy AC number to fam

+                  // if (!strncmp(fam,"PF",2)) strcut_(fam,'.'); // if PFAM identifier contains '.' cut it off

+                  // strcut_(ptr2); // cut after first word ...

+                  strcpy(name,ptr1); // ... and copy first word into name

+                  readCommentLine = '1'; /* MR1 */

+              }

+

+          //line contains sequence residues or SS information and does not belong to a >aa_ sequence

+          else if (!skip_sequence)

+              {

+                  if (v>=4) cout<<line<<"\n"; //DEBUG

+                  if (k==-1 && v)

+                      {

+                          cerr<<endl<<"WARNING: No sequence name preceding following line in "<<infile<<":\n\'"<<line<<"\'\n";

+                          continue;

+                      }

+

+                  h=0; //counts characters in current line

+

+                  // Check whether all characters are correct; store into cur_seq

+                  if (keep[k] || (k == kfirst) ) // normal line containing residues /* MR1 */

+                      {

+                          while (h<LINELEN && line[h]>'\0' && l</*MAXCOL*/par.maxColCnt-1)

+                              {

+                                  if (aa2i(line[h])>=0) // ignore white-space characters ' ', \t and \n (aa2i()==-1)

+                                      {cur_seq[l]=line[h]; l++;}

+                                  else if (aa2i(line[h])==-2 && v)

+                                      cerr<<endl<<"WARNING: invalid symbol \'"<<line[h]<<"\' at pos. "<<h<<" in line "<<linenr<<" of "<<infile<<"\n";

+                                  h++;

+                              }

+                      }

+                  else if (k==kss_dssp) // lines with dssp secondary structure states (. - H E C S T G B)

+                      {

+                          while (h<LINELEN && line[h]>'\0' && l</*MAXCOL*/par.maxColCnt-1)

+                              {

+                                  if (ss2i(line[h])>=0 && ss2i(line[h])<=7)

+                                      {cur_seq[l]=ss2ss(line[h]); l++;}

+                                  else if (v)

+                                      cerr<<endl<<"WARNING: invalid symbol \'"<<line[h]<<"\' at pos. "<<h<<" in line "<<linenr<<" of "<<infile<<"\n";

+                                  h++;

+                              }

+                      }

+                  else if (k==ksa_dssp) // lines with dssp solvent accessibility states (. - ???)

+                      {

+                          while (h<LINELEN && line[h]>'\0' && l</*MAXCOL*/par.maxColCnt-1)

+                              {

+                                  if (sa2i(line[h])>=0)

+                                      cur_seq[l++]=line[h];

+                                  else if (v)

+                                      cerr<<endl<<"WARNING: invalid symbol \'"<<line[h]<<"\' at pos. "<<h<<" in line "<<linenr<<" of "<<infile<<"\n";

+                                  h++;

+                              }

+                      }

+                  else if (k==kss_pred) // lines with predicted secondary structure (. - H E C)

+                      {

+                          while (h<LINELEN && line[h]>'\0' && l</*MAXCOL*/par.maxColCnt-1)

+                              {

+                                  if (ss2i(line[h])>=0 && ss2i(line[h])<=3)

+                                      {cur_seq[l]=ss2ss(line[h]); l++;}

+                                  else if (v)

+                                      cerr<<endl<<"WARNING: invalid symbol \'"<<line[h]<<"\' at pos. "<<h<<" in line "<<linenr<<" of "<<infile<<"\n";

+                                  h++;

+                              }

+                      }

+                  else if (k==kss_conf) // lines with confidence values should contain only 0-9, '-', or '.'

+                      {

+                          while (h<LINELEN && line[h]>'\0' && l</*MAXCOL*/par.maxColCnt-1)

+                              {

+                                  if (line[h]=='-' || line[h]=='.' || (line[h]>='0' && line[h]<='9'))

+                                      {cur_seq[l]=line[h]; l++;}

+                                  else if (v)

+                                      cerr<<endl<<"WARNING: invalid symbol \'"<<line[h]<<"\' at pos. "<<l<<" in line "<<linenr<<" of "<<infile<<"\n";

+                                  h++;

+                              }

+                      }

+                  else if (display[k]) // other lines such as >sa_pred etc

+                      {

+                          while (h<LINELEN && line[h]>'\0' && l</*MAXCOL*/par.maxColCnt-1)

+                              {

+                                  if (line[h]=='-' || line[h]=='.' || (line[h]>='0' && line[h]<='9') || (line[h]>='A' && line[h]<='B'))

+                                      {cur_seq[l]=line[h]; l++;}

+                                  else if (v)

+                                      cerr<<endl<<"WARNING: invalid symbol \'"<<line[h]<<"\' at pos. "<<l<<" in line "<<linenr<<" of "<<infile<<"\n";

+                                  h++;

+                              }

+                      }

+                  if (v && l>=/*MAXCOL*/par.maxColCnt-1) 

+                      {

+                          cerr<<endl<<"WARNING: maximum number of residues "<</*MAXCOL*/par.maxColCnt-2<<" exceded in sequence "<<sname[k]<<"\n";

+                          skip_sequence=1;

+                      }

+                  cur_seq[l]='\0'; //Ensure that cur_seq ends with a '\0' character

+              } //end else

+

+      }

+  /////////////////////////////////////////////////////////////////////////

+

+

+  if (k>=0) //if at least one sequence was read in

+      {

+          seq[k]=new(char[strlen(cur_seq)+2]);

+          if (!seq[k]) MemoryError("array for input sequences");

+          X[k]=new(char[strlen(cur_seq)+2]);

+          if (!X[k]) MemoryError("array for input sequences");

+          I[k]=new(short unsigned int[strlen(cur_seq)+2]);

+          if (!I[k]) MemoryError("array for input sequences");

+          strcpy(seq[k],cur_seq);

+      }

+  else

+      {cerr<<endl<<"Error: no sequences found in file "<<infile<<"\n"; throw 1;}

+

+  N_in = k+1;

+

+  // Set name, longname, fam

+  if (!*name) // longname, name and family were not set by '#...' line yet -> extract from first sequence

+      {

+          char* ptr;

+          // strtr(sname[kfirst],"~"," "); // 'transpose': replaces the tilde with a blanc everywhere in sname[kfirst]

+          strncpy(longname,sname[kfirst],DESCLEN-1); // longname is name of first sequence

+          longname[DESCLEN-1]='\0';

+          strncpy(name,sname[kfirst],NAMELEN-1); // Shortname is first word of longname...

+          name[NAMELEN-1]='\0';

+          ptr = strcut(name); // ...until first white-space character

+          if (ptr && islower(ptr[0]) && ptr[1]=='.' && isdigit(ptr[2])) //Scop family code present as second word?

+              {

+                  lwrstr(name); // Transform upper case to lower case

+                  strcut(ptr); // Non-white-space characters until next white-space character..

+                  strcpy(fam,ptr); // ...are the SCOP familiy code

+              }

+          else if (name[0]=='P' && name[1]=='F' && isdigit(name[2]) && isdigit(name[3]) ) //Pfam code

+              {

+                  strcpy(fam,name); // set family name = Pfam code

+              }

+      }

+  

+  

+  

+  delete[] cur_seq; cur_seq = NULL;

+  

+  // Checking for warning messages

+  if (v==0) return;

+  if (v>=2) cout<<"Read "<<infile<<" with "<<N_in<<" sequences\n";

+  if (v>=3) cout<<"Query sequence for alignment has number "<<kfirst<<" (0 is first)\n";

+  return;

+}

+

+/*

+ * At this point GetSeqsFromHMM() slots in, however,

+ * only needed in hhbliys.C, so will skip it for moment, MR1

+ */

+

+

+/////////////////////////////////////////////////////////////////////////////

+/**

+ * @brief  Convert ASCII in seq[k][l] to int (0-20) in X[k][i],

+ *  throw out all insert states, record their number in I[k][i]

+ *  and store sequences to be displayed in seq[k] */

+/////////////////////////////////////////////////////////////////////////////

+void 

+Alignment::Compress(const char infile[])

+{

+    int i; // Index for match state (first=1)

+    int l; // Postion in alignment incl. gaps (first=1)

+    int k; // Index for sequences (first=0)

+    int a; // amino acid index

+    char c;

+    int unequal_lengths=0; /* k: seq k doesn't have same number

+                              of match states as seq 0 => WARNING */

+    /* points to next character in seq[k] to be written */

+    /*static short unsigned int h[MAXSEQ];*/

+    /*short*/ unsigned int *h = NULL; /* short may lead to overflow for long alignments, FS, r235 -> r236 */

+

+    h = new(/*short*/ unsigned int[N_in+2]); /* short -> overflow, FS, r235 -> r236 */

+    float *percent_gaps = NULL; /* FS, 2010-Nov */

+    char *match_state = NULL;  /* FS, 2010-Nov */

+

+

+    // Initialize 

+    for (k=0;k<N_in; k++) 

+        {I[k][0]=0;}

+

+    if (v>=3)

+        {

+            if (par.M==1)

+                cout<<"Using match state assignment by capital letters (a2m format)\n";

+            else if (par.M==2) cout<<"Using percentage-rule match state assignment\n";

+            else if (par.M==3) cout<<"Using residues of first sequence as match states\n";

+        }

+

+    // Create matrices X and I with amino acids represented by integer numbers

+    switch(par.M)

+        {

+

+            /////////////////////////////////////////////////////////////////////////

+            /* a2m/a3m format: match states capital case,

+               inserts lower case, delete states '-', inserted gaps '.'

+               The confidence values for ss prediction are interpreted as follows:

+               0-9:match states(!) '-' :match state '.':insert */

+        case 1:

+        default:

+

+            // Warn if alignment is ment to be -M first or -M NN instead of A2M/A3M

+            if (v>=2 && strchr(seq[kfirst],'-') ) // Seed/query sequence contains a gap ...

+                {

+                    for (k=1; k<N_in; k++)

+                        if (strpbrk(seq[k],"abcdefghiklmnpqrstuvwxyz.")) break;

+                    if (k==N_in) // ... but alignment contains no lower case residue

+                        printf("WARNING: input alignment %s looks like aligned FASTA instead of A2M/A3M format. Consider using '-M first' or '-M 50'\n",infile);

+                }

+

+            // Remove '.' characters from seq[k]

+            for(k=0; k<N_in; k++)

+                {

+                    char* ptrS=seq[k]; // pointer to source: character in seq[k]

+                    char* ptrD=seq[k]; // pointer to destination: seq[k]

+                    while(1) // omit '.' symbols

+                        {

+                            if (*ptrS!='.') {*ptrD=*ptrS; ptrD++;} //leave out '.' symbols

+                            if (!*ptrS) break;

+                            ptrS++;

+                        }

+                }

+            L=/*MAXRES*/par.maxResLen-2; // needed because L=imin(L,i)

+            for (k=0; k<N_in; k++)

+                {

+                    i=1; l=1; // start at i=1, not i=0!

+                    if (keep[k]) //skip >ss_dssp, >ss_pred, >ss_conf, >aa_... sequences

+                        {

+                            while((c=seq[k][l++])) // assign residue to c at same time

+                                {

+                                    if (c>='a' && c<='z') I[k][i-1]++;//insert state = lower case character

+                                    else if (c!='.') //match state = upper case character

+                                        {

+                                            X[k][i]=aa2i(c);

+                                            I[k][i]=0;

+                                            i++;

+                                        }

+                                }

+                        }

+                    else if (k==kss_dssp || k==kss_pred) // does alignment contain sequence of secondary structure states?

+                        {

+                            while((c=seq[k][l++])) // assign residue to c at same time

+                                if (c!='.' && !(c>='a' && c<='z')) X[k][i++]=ss2i(c); //match state = upper case character

+                        }

+                    else if (k==ksa_dssp) // does alignment contain sequence of prediction confidence values?

+                        {

+                            while((c=seq[k][l++])) // assign residue to c at same time

+                                if (c!='.' && !(c>='a' && c<='z')) X[k][i++]=sa2i(c); //match state = upper case character

+                        }

+                    else if (k==kss_conf) // does alignment contain sequence of prediction confidence values?

+                        {

+                            while((c=seq[k][l++])) // assign residue to c at same time

+                                if (c!='.') X[k][i++]=cf2i(c); //match state = 0-9 or '-'

+                        }

+                    else if (k==kfirst)        // does alignment contain sequence of prediction confidence values?

+                        {

+                            while((c=seq[k][l++]))  // assign residue to c at same time