binaries/src/clustalo/src/hhalign/hhutil-C.h

   1 /* -*- mode: c; tab-width: 4; c-basic-offset: 4; indent-tabs-mode: nil -*- */
   2
   3 /*********************************************************************
   4  * Clustal Omega - Multiple sequence alignment
   5  *
   6  * Copyright (C) 2010 University College Dublin
   7  *
   8  * Clustal-Omega is free software; you can redistribute it and/or
   9  * modify it under the terms of the GNU General Public License as
  10  * published by the Free Software Foundation; either version 2 of the
  11  * License, or (at your option) any later version.
  12  *
  13  * This file is part of Clustal-Omega.
  14  *
  15  ********************************************************************/
  16
  17 /*
  18  * RCS $Id: hhutil-C.h 246 2011-06-15 12:41:04Z fabian $
  19  */
  20
  21 //////////////////////////////////////////////////////////////////////////////
  22 // Transform a character to lower case and '.' to '-' and vice versa
  23 //////////////////////////////////////////////////////////////////////////////
  24
  25
  26 inline char
  27 MatchChr(char c)  {return ((c>='a' && c<='z')? c-'a'+'A' : (c=='.'? '-':c) );}
  28
  29 inline char
  30 InsertChr(char c) {return ((c>='A' && c<='Z')? c+'a'-'A' : ((c>='0' && c<='9') || c=='-')? '.':c );}
  31
  32 inline int
  33 WordChr(char c) {return (int)((c>='A' && c<='Z') || (c>='a' && c<='z'));}
  34
  35
  36 //////////////////////////////////////////////////////////////////////////////
  37 /**
  38  * @brief Transforms the one-letter amino acid code into an integer between 0 and 22
  39  */
  40 inline char
  41 aa2i(char c)
  42 {
  43   //A  R  N  D  C  Q  E  G  H  I  L  K  M  F  P  S  T  W  Y  V
  44   if (c>='a' && c<='z') c+='A'-'a';
  45   switch (c)
  46     {
  47     case 'A': return 0;
  48     case 'R': return 1;
  49     case 'N': return 2;
  50     case 'D': return 3;
  51     case 'C': return 4;
  52     case 'Q': return 5;
  53     case 'E': return 6;
  54     case 'G': return 7;
  55     case 'H': return 8;
  56     case 'I': return 9;
  57     case 'L': return 10;
  58     case 'K': return 11;
  59     case 'M': return 12;
  60     case 'F': return 13;
  61     case 'P': return 14;
  62     case 'S': return 15;
  63     case 'T': return 16;
  64     case 'W': return 17;
  65     case 'Y': return 18;
  66     case 'V': return 19;
  67     case 'X': return ANY;
  68     case 'J': return ANY;
  69     case 'O': return ANY;
  70     case 'U': return 4;  //Selenocystein -> Cystein
  71     case 'B': return 3;  //D (or N)
  72     case 'Z': return 6;  //E (or Q)
  73     case '-': return GAP;
  74     case '.': return GAP;
  75     case '_': return GAP;
  76     }
  77   if (c>=0 && c<=32) return -1; // white space and control characters
  78   return -2;
  79 }
  80
  81 ///////////////////////////////////////////////////////////////////////////////
  82 /**
  83  * @brief Transforms integers between 0 and 22 into the one-letter amino acid code
  84  */
  85 inline char
  86 i2aa(char c)
  87 {
  88   //A  R  N  D  C  Q  E  G  H  I  L  K  M  F  P  S  T  W  Y  V
  89   switch (c)
  90     {
  91     case 0: return 'A';
  92     case 1: return 'R';
  93     case 2: return 'N';
  94     case 3: return 'D';
  95     case 4: return 'C';
  96     case 5: return 'Q';
  97     case 6: return 'E';
  98     case 7: return 'G';
  99     case 8: return 'H';
 100     case 9: return 'I';
 101     case 10: return 'L';
 102     case 11: return 'K';
 103     case 12: return 'M';
 104     case 13: return 'F';
 105     case 14: return 'P';
 106     case 15: return 'S';
 107     case 16: return 'T';
 108     case 17: return 'W';
 109     case 18: return 'Y';
 110     case 19: return 'V';
 111     case ANY: return 'X';
 112     case GAP: return '-';
 113     case ENDGAP: return '-';
 114     }
 115   return '?';
 116 }
 117
 118 //////////////////////////////////////////////////////////////////////////////
 119 /**
 120  * @brief Transforms the dssp/psipred secondary structure code into an integer number
 121  */
 122 inline char
 123 ss2i(char c)
 124 {
 125   //- H E C S T G B
 126   if (c>='a' && c<='z') c+='A'-'a';
 127   switch (c)
 128     {
 129     case '.': return 0;
 130     case '-': return 0;
 131     case 'X': return 0;
 132     case 'H': return 1;
 133     case 'E': return 2;
 134     case 'C': return 3;
 135     case '~': return 3;
 136     case 'S': return 4;
 137     case 'T': return 5;
 138     case 'G': return 6;
 139     case 'B': return 7;
 140     case 'I': return 3;
 141     case ' ': return -1;
 142     case '\t': return -1;
 143     case '\n': return -1;
 144     }
 145   return -2;
 146 }
 147
 148 //////////////////////////////////////////////////////////////////////////////
 149 /**
 150  * @brief Transforms integers between 0 and 8 into the dssp/psipred secondary structure code
 151  */
 152 inline char
 153 i2ss(int c)
 154 {
 155   //- H E C S T G B
 156   switch (c)
 157     {
 158     case 0: return '-';
 159     case 1: return 'H';
 160     case 2: return 'E';
 161     case 3: return 'C';
 162     case 4: return 'S';
 163     case 5: return 'T';
 164     case 6: return 'G';
 165     case 7: return 'B';
 166     case 8: return 'I';
 167     }
 168   return '?';
 169 }
 170
 171
 172 //////////////////////////////////////////////////////////////////////////////
 173 /**
 174  * @brief Transforms the solvend accessiblity code into an integer number
 175  */
 176 inline char
 177 sa2i(char c)
 178 {
 179   //- A B C D E
 180   if (c>='a' && c<='z') c+='A'-'a';
 181   switch (c)
 182     {
 183     case '.': return 0;
 184     case '-': return 0;
 185     case 'A': return 1;
 186     case 'B': return 2;
 187     case 'C': return 3;
 188     case 'D': return 4;
 189     case 'E': return 5;
 190     case 'F': return 6;
 191     case ' ': return -1;
 192     case '\t': return -1;
 193     case '\n': return -1;
 194     }
 195   return -2;
 196 }
 197
 198 //////////////////////////////////////////////////////////////////////////////
 199 /**
 200  * @brief Transforms integers between 0 and 5 into the solvent accessibility code
 201  */
 202 inline char i2sa(int c)
 203 {
 204   //- H E C S T G B
 205   switch (c)
 206     {
 207     case 0: return '-';
 208     case 1: return 'A';
 209     case 2: return 'B';
 210     case 3: return 'C';
 211     case 4: return 'D';
 212     case 5: return 'E';
 213     case 6: return 'F';
 214     }
 215   return '?';
 216 }
 217
 218
 219 //////////////////////////////////////////////////////////////////////////////
 220 /**
 221  * @brief Transforms alternative secondary structure symbols into symbols
 222  */
 223 inline char
 224 ss2ss(char c)
 225 {
 226   //- H E C S T G B
 227   switch (c)
 228     {
 229     case '~': return 'C';
 230     case 'I': return 'C';
 231     case 'i': return 'c';
 232     case 'H':
 233     case 'E':
 234     case 'C':
 235     case 'S':
 236     case 'T':
 237     case 'G':
 238     case 'B':
 239     case 'h':
 240     case 'e':
 241     case 'c':
 242     case 's':
 243     case 't':
 244     case 'g':
 245     case 'b':
 246     case '.':
 247       return c;
 248     }
 249   return '-';
 250 }
 251
 252 //////////////////////////////////////////////////////////////////////////////
 253 /**
 254  * @brief Transforms confidence values of psipred into internal code
 255  */
 256 inline char
 257 cf2i(char c)
 258 {
 259   switch (c)
 260     {
 261     case '-': return 0;
 262     case '.': return 0;
 263     case '0': return 1;
 264     case '1': return 2;
 265     case '2': return 3;
 266     case '3': return 4;
 267     case '4': return 5;
 268     case '5': return 6;
 269     case '6': return 7;
 270     case '7': return 8;
 271     case '8': return 9;
 272     case '9': return 10;
 273     }
 274   return 0;
 275 }
 276
 277 //////////////////////////////////////////////////////////////////////////////
 278 /**
 279  * @brief Transforms internal representation of psipred confidence values into printable chars
 280  */
 281 inline char
 282 i2cf(char c)
 283 {
 284   switch (c)
 285     {
 286     case 0: return '-';
 287     case 1: return '0';
 288     case 2: return '1';
 289     case 3: return '2';
 290     case 4: return '3';
 291     case 5: return '4';
 292     case 6: return '5';
 293     case 7: return '6';
 294     case 8: return '7';
 295     case 9: return '8';
 296     case 10: return '9';
 297     }
 298   return '-';
 299 }
 300
 301
 302 //////////////////////////////////////////////////////////////////////////////
 303 /**
 304  * @brief Fast lookup of log2(1+2^(-x)) for x>=0 (precision < 0.35%)
 305  */
 306 inline float
 307 fast_addscore(float x)
 308 {
 309   static float val[2001];         // val[i]=log2(1+2^(-x))
 310   static char initialized;
 311   if (x>20) return 0.0;
 312   if (x<0)
 313     {
 314       fprintf(stderr,"Error in function fast_addscore: argument %g is negative\n",x);
 315       exit(7);
 316     }
 317   if (!initialized)   //First fill in the log2-vector
 318     {
 319       for (int i=0; i<=2000; i++) val[i]=log2(1.0+pow(2,-0.01*(i+0.5)));
 320       initialized=1;
 321     }
 322   return val[(int)(100.0*x)];
 323 }
 324
 325
 326
 327 //////////////////////////////////////////////////////////////////////////////
 328 /**
 329  * @brief Little utilities for output
 330  */
 331 inline void
 332 fout(FILE* outf, int d)
 333 {
 334   if (d>=99999) fprintf(outf,"*\t"); else fprintf(outf,"%i\t",d);
 335   return;
 336 }
 337
 338 //////////////////////////////////////////////////////////////////////////////
 339 /**
 340  * @brief Errors
 341  */
 342 int
 343 FormatError(const char infile[], const char details[]="")
 344 {
 345   cerr<<"Error in "<</*par.argv[0],FS*/__FILE__<<": wrong format while reading file \'"<<infile<<". "<<details<<"\n";
 346   exit(1);
 347 }
 348
 349 int
 350 OpenFileError(const char outfile[])
 351 {
 352   cerr<<endl<<"Error in "<</*par.argv[0],FS*/__FILE__<<": could not open file \'"<<outfile<<"\'\n";
 353   exit(2);
 354 }
 355
 356 int
 357 MemoryError(const char arrayname[])
 358 {
 359   cerr<<"Error in "<</*par.argv[0],FS*/__FILE__<<": Memory overflow while creating \'"<<arrayname<<"\'. Please report this bug to developers\n";
 360   exit(3);
 361 }
 362
 363 int
 364 InternalError(const char errstr[])
 365 {
 366   cerr<<"Error in "<</*par.argv[0],FS*/__FILE__<<":  "<<errstr<<". Please report this bug to developers\n";
 367   exit(6);
 368 }
 369
 370
 371 //////////////////////////////////////////////////////////////////////////////
 372 /**
 373  * @brief Takes family code (eg. a.1.2.3) and returns strings 'a', 'a.1', and 'a.1.2'
 374  */
 375 inline void
 376 ScopID(char cl[], char fold[], char sfam[], const char fam[])
 377 {
 378   char* ptr;
 379
 380   //get scop class ID
 381   strcpy(cl,fam);
 382   ptr = strchr(cl,'.');               //return adress of next '.' in name
 383   if(ptr) ptr[0]='\0';
 384
 385   //get scop fold ID
 386   strcpy(fold,fam);
 387   ptr = strchr(fold,'.');             //return adress of next '.' in name
 388   if(ptr) ptr = strchr(ptr+1,'.');    //return adress of next '.' in name
 389   if(ptr) ptr[0]='\0';
 390
 391   //get scop superfamily ID
 392   strcpy(sfam,fam);
 393   ptr = strchr(sfam,'.');            //return adress of next '.' in name
 394   if(ptr) ptr = strchr(ptr+1,'.');   //return adress of next '.' in name
 395   if(ptr) ptr = strchr(ptr+1,'.');   //return adress of next '.' in name
 396   if(ptr) ptr[0]='\0';
 397   return;
 398 }
 399
 400 //////////////////////////////////////////////////////////////////////////////
 401 /**
 402  * @brief Read up to n lines of outfile and write to screen (STDERR)
 403  */
 404 void
 405 WriteToScreen(char* outfile, int n)
 406 {
 407   char line[LINELEN]="";
 408   ifstream outf;
 409   outf.open(outfile, ios::in);
 410   if (!outf) {OpenFileError(outfile);}
 411   cout<<"\n";
 412   for(; n>0 && outf.getline(line,LINELEN); n--) cout<<line<<"\n";
 413   outf.close();
 414   cout<<"\n";
 415 }
 416
 417 inline void
 418 WriteToScreen(char* outfile) {WriteToScreen(outfile,INT_MAX);}
 419
 420
 421
 422 /////////////////////////////////////////////////////////////////////////////////////
 423 /**
 424  * @brief Read .hhdefaults file into array argv_conf (beginning at argv_conf[1])
 425  */
 426 void
 427 ReadDefaultsFile(int& argc_conf, char** argv_conf)
 428 {
 429   char line[LINELEN]="";
 430   char filename[NAMELEN];
 431   char* c_first;   //pointer to first character of argument string
 432   char* c;         //pointer to scan line read in for end of argument
 433   //  ifstream configf;
 434   FILE* configf=NULL;
 435   argc_conf=1;     //counts number of arguments read in
 436
 437   // Open config file
 438   strcpy(filename,"./.hhdefaults");
 439   configf = fopen(filename,"r");
 440   if (!configf && getenv("HOME"))
 441     {
 442       strcpy(filename,getenv("HOME"));
 443       strcat(filename,"/.hhdefaults");
 444       configf = fopen(filename,"r");
 445       if (!configf)
 446         {
 447           if (v>=3) cerr<<"Warning: could not find ./.hhdefaults or "<<filename<<"\n";
 448           return;
 449         }
 450     }
 451   else if (!configf) return; // only webserver has no home directory => need no warning
 452
 453   // Scan file until line 'program_nameANYTHING'
 454   while (fgets(line,LINELEN,configf))
 455     if (!strncmp(line,program_name,6)) break;
 456   // Found line 'program_nameANYTHING'?
 457   if (!strncmp(line,program_name,6))
 458     {
 459       // Read in options until end-of-file or empty line
 460         while (fgets(line,LINELEN,configf) && strcmp(line,"\n"))
 461         {
 462           // Analyze line
 463           c=line;
 464           do
 465             {
 466               // Find next word
 467               while (*c==' ' || *c=='\t') c++; //Advance until next non-white space
 468               if (*c=='\0' || *c=='\n' || *c=='#') break;  //Is next word empty string?
 469               c_first=c;
 470               while (*c!=' ' && *c!='\t'  && *c!='#' && *c!='\0' && *c!='\n' ) c++; //Advance until next white space or '#'
 471               if (*c=='\0' || *c=='\n' || *c=='#')         //Is end of line reached?
 472                 {
 473                   *c='\0';
 474                   argv_conf[argc_conf]=new(char[strlen(c_first)+1]);
 475                   strcpy(argv_conf[argc_conf++],c_first);
 476                   break;
 477                 }
 478               *c='\0';
 479               argv_conf[argc_conf]=new(char[strlen(c_first)+1]);
 480               strcpy(argv_conf[argc_conf++],c_first);
 481               printf("Argument: %s\n",c_first);
 482               c++;
 483             } while (1);
 484         } //end read line
 485      if (v>=3)
 486         {
 487           cout<<"Arguments read in from .hhdefaults:";
 488           for (int argc=1; argc<argc_conf; argc++) cout<<(argv_conf[argc][0]=='-'? " ":"")<<argv_conf[argc]<<" ";
 489           cout<<"\n";
 490         }
 491      else if (v>=3) cout<<"Read in "<<argc_conf<<" default arguments for "<<program_name<<" from "<<filename<<"\n";
 492      }
 493   else //found no line 'program_name   anything"
 494     {
 495       if (v>=3) cerr<<endl<<"Warning: no default options for \'"<<program_name<<"\' found in "<<filename<<"\n";
 496       return; //no line 'program_name   anything' found
 497     }
 498 //   configf.close();
 499   fclose(configf);
 500 }
 501
 502
 503 /////////////////////////////////////////////////////////////////////////////////////
 504 /**
 505  * @brief Set default parameter values
 506  */
 507 void
 508 SetDefaults()
 509 {
 510
 511   par.append=0;                // overwrite output file
 512   par.outformat=0;             // 0: hhr  1: FASTA  2:A2M   3:A3M
 513   par.p=20.0f;                 // minimum threshold for inclusion in hit list and alignment listing
 514   par.E=1e6f;                  // maximum threshold for inclusion in hit list and alignment listing
 515   par.b=10;                    // min number of alignments
 516   par.B=500;                   // max number of alignments
 517   par.z=10;                    // min number of lines in hit list
 518   par.Z=500;                   // max number of lines in hit list
 519   par.e=1e-3f;                 // maximum E-value for inclusion in output alignment, output HMM, and PSI-BLAST checkpoint model
 520   par.showcons=1;              // show consensus sequence
 521   par.showdssp=1;              // show predicted secondary structure
 522   par.showpred=1;              // show dssp secondary structure
 523   par.cons=0;                  // show first non-SS sequence as main representative sequence (not consensus)
 524   par.nseqdis=1;               // maximum number of query sequences for output alignment
 525   par.mark=0;                  // 1: only marked sequences (or first) get displayed; 0: most divergent ones get displayed
 526   par.aliwidth=80;             // number of characters per line in output alignments for HMM search
 527   par.max_seqid=90;            // default for maximum sequence identity threshold
 528   par.qid=0;                   // default for minimum sequence identity with query
 529   par.qsc=-20.0f;              // default for minimum score per column with query
 530   par.coverage=0;              // default for minimum coverage threshold
 531   par.Ndiff=100;               // pick Ndiff most different sequences from alignment
 532   par.coverage_core=80;        // Minimum coverage for sequences in core alignment
 533   par.qsc_core=0.3f;           // Minimum score per column of core sequence with query
 534   par.coresc=-20.0f;           // Minimum score per column with core alignment (HMM)
 535
 536   par.M=1;                     // match state assignment is by A2M/A3M
 537   par.Mgaps=50;                // Above this percentage of gaps, columns are assigned to insert states (for par.M=2)
 538   par.calibrate=0;             // default: no calibration
 539   par.calm=0;                  // derive P-values from: 0:query calibration  1:template calibration  2:both
 540   par.mode=0;                  //
 541
 542   par.wg=0;                    // 0: use local sequence weights   1: use local ones
 543
 544   par.matrix=0;                // Subst.matrix 0: Gonnet, 1: HSDM, 2: BLOSUM50 3: BLOSUM62
 545   par.pcm=2;                   // pseudocount mode: default=divergence-dependent (but not column-specific)
 546 #if 1 /* Nelder-Meade on Baliscore */
 547   par.pca=1.712190f;                // default values for substitution matrix pseudocounts
 548   par.pcb=1.039640f;                // significant reduction of pcs by Neff_M starts around Neff_M-1=pcb
 549   par.pcc=0.878067f;                // pcs are reduced prop. to 1/Neff^pcc
 550   par.pcw=0.0f;                // wc>0 weighs columns according to their intra-clomun similarity
 551
 552   par.gapb=1.405220;                // default values for transition pseudocounts
 553   par.gapd=1.316760;               // gap open penalty pseudocount; 0.25 corresponds to 7.1*gapf bits
 554   par.gape=1.793780;                // gap extension penalty pseudocount
 555   par.gapf=1.034710;                // factor for increasing gap open penalty for deletes
 556   par.gapg=0.894772;                // factor for increasing gap open penalty for inserts
 557   par.gaph=0.544072;                // factor for increasing gap extension penalty for deletes
 558   par.gapi=0.862559;                // factor for increasing gap extension penalty for inserts
 559 #else  /* Soeding's default*/
 560   par.pca=1.0f;                // default values for substitution matrix pseudocounts
 561   par.pcb=1.5f;                // significant reduction of pcs by Neff_M starts around Neff_M-1=pcb
 562   par.pcc=1.0f;                // pcs are reduced prop. to 1/Neff^pcc
 563   par.pcw=0.0f;                // wc>0 weighs columns according to their intra-clomun similarity
 564
 565   par.gapb=1.0;                // default values for transition pseudocounts
 566   par.gapd=0.15;               // gap open penalty pseudocount; 0.25 corresponds to 7.1*gapf bits
 567   par.gape=1.0;                // gap extension penalty pseudocount
 568   par.gapf=0.6;                // factor for increasing gap open penalty for deletes
 569   par.gapg=0.6;                // factor for increasing gap open penalty for inserts
 570   par.gaph=0.6;                // factor for increasing gap extension penalty for deletes
 571   par.gapi=0.6;                // factor for increasing gap extension penalty for inserts
 572 #endif
 573
 574 #if 0
 575   /* Viterbi parameters optimised on Sabre (R228), FS, r228 -> r229 */
 576   par.pcaV=1.245150f;                // default values for substitution matrix pseudocounts
 577   par.pcbV=1.682110f;                // significant reduction of pcs by Neff_M starts around Neff_M-1=pcb
 578   par.pccV=1.483840f;                // pcs are reduced prop. to 1/Neff^pcc
 579   par.pcwV=0.0f;                // wc>0 weighs columns according to their intra-clomun similarity
 580
 581   par.gapbV=0.818625;                // default values for transition pseudocounts
 582   par.gapdV=0.666110;               // gap open penalty pseudocount; 0.25 corresponds to 7.1*gapf bits
 583   par.gapeV=1.028050;                // gap extension penalty pseudocount
 584   par.gapfV=0.710760;                // factor for increasing gap open penalty for deletes
 585   par.gapgV=1.649800;                // factor for increasing gap open penalty for inserts
 586   par.gaphV=0.470604;                // factor for increasing gap extension penalty for deletes
 587   par.gapiV=0.829479;                // factor for increasing gap extension penalty for inserts
 588 #elif 1
 589   /* Viterbi parameters optimised on Balibase, r244 -> r245 */
 590   par.pcaV=1.333860f;                // default values for substitution matrix pseudocounts
 591   par.pcbV=1.934480f;                // significant reduction of pcs by Neff_M starts around Neff_M-1=pcb
 592   par.pccV=1.655610f;                // pcs are reduced prop. to 1/Neff^pcc
 593   par.pcwV=0.0f;                // wc>0 weighs columns according to their intra-clomun similarity
 594
 595   par.gapbV=0.334525;                // default values for transition pseudocounts
 596   par.gapdV=0.074534;               // gap open penalty pseudocount; 0.25 corresponds to 7.1*gapf bits
 597   par.gapeV=0.320336;                // gap extension penalty pseudocount
 598   par.gapfV=0.151634;                // factor for increasing gap open penalty for deletes
 599   par.gapgV=0.641516;                // factor for increasing gap open penalty for inserts
 600   par.gaphV=0.266434;                // factor for increasing gap extension penalty for deletes
 601   par.gapiV=0.598414;                // factor for increasing gap extension penalty for inserts
 602 #else /* Soeding default*/
 603   par.pcaV=1.0f;                // default values for substitution matrix pseudocounts
 604   par.pcbV=1.5f;                // significant reduction of pcs by Neff_M starts around Neff_M-1=pcb
 605   par.pccV=1.0f;                // pcs are reduced prop. to 1/Neff^pcc
 606   par.pcwV=0.0f;                // wc>0 weighs columns according to their intra-clomun similarity
 607
 608   par.gapbV=1.0;                // default values for transition pseudocounts
 609   par.gapdV=0.15;               // gap open penalty pseudocount; 0.25 corresponds to 7.1*gapf bits
 610   par.gapeV=1.0;                // gap extension penalty pseudocount
 611   par.gapfV=0.6;                // factor for increasing gap open penalty for deletes
 612   par.gapgV=0.6;                // factor for increasing gap open penalty for inserts
 613   par.gaphV=0.6;                // factor for increasing gap extension penalty for deletes
 614   par.gapiV=0.6;                // factor for increasing gap extension penalty for inserts
 615 #endif
 616
 617   par.ssm=2;                   // ss scoring mode: 0:no ss score  1:score after alignment  2:score during alignment
 618   par.ssw=0.11f;               // weight of ss scoring
 619   par.ssa=1.0f;                // weight of ss evolution matrix
 620   par.shift=-0.01f;            // Shift match score up
 621   par.mact=0.3001f;            // Score threshold for MAC alignment in local mode (set to 0.5001 to track user modification)
 622   par.corr=0.1f;               // Weight of correlations of scores for |i-j|<=4
 623   par.wstruc=1.0f;             // Weight of structure scores
 624
 625   par.egq=0.0f;                // no charge for end gaps as default
 626   par.egt=0.0f;                // no charge for end gaps as default
 627
 628   par.trans=0;                 // no transitive scoring as default
 629   par.Emax_trans=100.0f;       // use intermediate HMMs with E-values up to 100 between query and database HMM
 630   par.Emax_trans=100.0f;       // use intermediate HMMs with E-values up to 100 between query and database HMM
 631   par.wtrans=1.0f;             // Ztot[k] = Zq[k] + wtrans * (Zforward[k]+Zreverse[k])
 632   par.ssgap=0;                 // 1: add secondary structure-dependent gap penalties  0:off
 633   par.ssgapd=1.0f;             // secondary structure-dependent gap-opening penalty (per residue)
 634   par.ssgape=0.0f;             // secondary structure-dependent gap-extension penalty (per residue)
 635   par.ssgapi=4;                // max. number of inside-integer(ii); gap-open-penalty= -ii*ssgapd
 636
 637   par.loc=1;                   // local vs. global alignment as default
 638   par.altali=2;                // find up to two (possibly overlapping) subalignments
 639   par.forward=0;               // 0: Viterbi algorithm; 1: Viterbi+stochastic sampling; 3:Maximum Accuracy (MAC) algorithm
 640   par.realign=1;               // realign with MAC algorithm
 641
 642   par.repmode=0;               // repeats score independently of one another
 643   par.columnscore=1;           // Default column score is 1: null model pnul = 1/2 * (q_av(a)+p_av(a))
 644   par.min_overlap=0;           // automatic minimum overlap used
 645   par.opt=0;                   // Default = optimization mode off
 646   par.readdefaultsfile=0;      // Default = do not read a defaults file ./.hhdefaults or HOME/.hhdefaults
 647   par.maxdbstrlen=200;         // maximum length of database string to be printed in 'Command' line of hhr file
 648   par.mode=0;
 649   par.idummy=par.jdummy=0;     //
 650
 651   par.notags=1;                // neutralize His-tags, FLAG-tags, C-myc-tags
 652
 653   // Initialize strings
 654   strcpy(par.infile,"stdin");
 655   strcpy(par.outfile,"");
 656   strcpy(par. pairwisealisfile,"");
 657   strcpy(par.buffer,"buffer.txt");
 658   strcpy(par.scorefile,"");
 659   strcpy(par.wfile,"");
 660   strcpy(par.alnfile,"");
 661   strcpy(par.hhmfile,"");
 662   strcpy(par.psifile,"");
 663   par.exclstr=NULL;
 664
 665 #if 0  /* read parameter file from home-dir */
 666 #include "hhutil-C-help.h"
 667 #endif /* read parameter file from home-dir */
 668
 669    return;
 670 } /** this is the end of SetDefaults() **/
 671
 672 /*
 673  * EOF hhutil-C.h
 674  */