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 315 2016-12-15 17:18:30Z fabian $
  19  */
  20
  21 #define PARAMETERSFROMFILE 0
  22
  23 extern bool nucleomode;
  24
  25 //////////////////////////////////////////////////////////////////////////////
  26 // Transform a character to lower case and '.' to '-' and vice versa
  27 //////////////////////////////////////////////////////////////////////////////
  28
  29 inline char
  30 MatchChr(char c)  {return ((c>='a' && c<='z')? c-'a'+'A' : (c=='.'? '-':c) );}
  31
  32 inline char
  33 InsertChr(char c) {return ((c>='A' && c<='Z')? c+'a'-'A' : ((c>='0' && c<='9') || c=='-')? '.':c );}
  34
  35 inline int
  36 WordChr(char c) {return (int)((c>='A' && c<='Z') || (c>='a' && c<='z'));}
  37
  38
  39 //////////////////////////////////////////////////////////////////////////////
  40 /**
  41  * @brief Transforms the one-letter amino acid code into an integer between 0 and 22
  42  */
  43 inline char
  44 aa2i(char c)
  45 {
  46   if (c>='a' && c<='z') c+='A'-'a';
  47   if(!nucleomode)
  48   {
  49           // proteins!
  50           // A  R  N  D  C  Q  E  G  H  I  L  K  M  F  P  S  T  W  Y  V
  51           switch (c)
  52                 {
  53                 case 'A': return 0;
  54                 case 'R': return 1;
  55                 case 'N': return 2;
  56                 case 'D': return 3;
  57                 case 'C': return 4;
  58                 case 'Q': return 5;
  59                 case 'E': return 6;
  60                 case 'G': return 7;
  61                 case 'H': return 8;
  62                 case 'I': return 9;
  63                 case 'L': return 10;
  64                 case 'K': return 11;
  65                 case 'M': return 12;
  66                 case 'F': return 13;
  67                 case 'P': return 14;
  68                 case 'S': return 15;
  69                 case 'T': return 16;
  70                 case 'W': return 17;
  71                 case 'Y': return 18;
  72                 case 'V': return 19;
  73                 case 'X': return ANY;
  74                 case 'J': return ANY;
  75                 case 'O': return ANY;
  76                 case 'U': return 4;  //Selenocystein -> Cystein
  77                 case 'B': return 3;  //D (or N)
  78                 case 'Z': return 6;  //E (or Q)
  79                 case '-': return GAP;
  80                 case '.': return GAP;
  81                 case '_': return GAP;
  82                 }
  83   }
  84   else
  85   {
  86                 // nucleotides!
  87         switch(c)
  88         {
  89     case 'A': return 0;
  90     case 'C': return 1;
  91     case 'G': return 2;
  92     case 'T': return 3;
  93     case 'U': return 4;
  94     case 'N': return ANY;
  95     case '-': return GAP;
  96     case '.': return GAP;
  97     case '_': return GAP;
  98         default: return ANY;
  99         }
 100   }
 101   if (c>=0 && c<=32) return -1; // white space and control characters
 102   return -2;
 103 }
 104
 105 ///////////////////////////////////////////////////////////////////////////////
 106 /**
 107  * @brief Transforms integers between 0 and 22 into the one-letter amino acid code
 108  */
 109 inline char
 110 i2aa(char c)
 111 {
 112   //A  R  N  D  C  Q  E  G  H  I  L  K  M  F  P  S  T  W  Y  V
 113   if(!nucleomode)
 114   {
 115   switch (c)
 116     {
 117     case 0: return 'A';
 118     case 1: return 'R';
 119     case 2: return 'N';
 120     case 3: return 'D';
 121     case 4: return 'C';
 122     case 5: return 'Q';
 123     case 6: return 'E';
 124     case 7: return 'G';
 125     case 8: return 'H';
 126     case 9: return 'I';
 127     case 10: return 'L';
 128     case 11: return 'K';
 129     case 12: return 'M';
 130     case 13: return 'F';
 131     case 14: return 'P';
 132     case 15: return 'S';
 133     case 16: return 'T';
 134     case 17: return 'W';
 135     case 18: return 'Y';
 136     case 19: return 'V';
 137     case ANY: return 'X';
 138     case GAP: return '-';
 139     case ENDGAP: return '-';
 140     }
 141   }
 142   else
 143   {
 144         // DNA/RNA
 145         switch(c)
 146         {
 147                 case 0: return 'A';
 148                 case 1: return 'C';
 149                 case 2: return 'G';
 150                 case 3: return 'T';
 151                 case 4: return 'U';
 152                 case ANY: return 'N';
 153                 case GAP: return '-';
 154                 case ENDGAP: return '-';
 155         }
 156   }
 157   return '?';
 158 }
 159
 160 //////////////////////////////////////////////////////////////////////////////
 161 /**
 162  * @brief Transforms the dssp/psipred secondary structure code into an integer number
 163  */
 164 inline char
 165 ss2i(char c)
 166 {
 167   //- H E C S T G B
 168   if (c>='a' && c<='z') c+='A'-'a';
 169   switch (c)
 170     {
 171     case '.': return 0;
 172     case '-': return 0;
 173     case 'X': return 0;
 174     case 'H': return 1;
 175     case 'E': return 2;
 176     case 'C': return 3;
 177     case '~': return 3;
 178     case 'S': return 4;
 179     case 'T': return 5;
 180     case 'G': return 6;
 181     case 'B': return 7;
 182     case 'I': return 3;
 183     case ' ': return -1;
 184     case '\t': return -1;
 185     case '\n': return -1;
 186     }
 187   return -2;
 188 }
 189
 190 //////////////////////////////////////////////////////////////////////////////
 191 /**
 192  * @brief Transforms integers between 0 and 8 into the dssp/psipred secondary structure code
 193  */
 194 inline char
 195 i2ss(int c)
 196 {
 197   //- H E C S T G B
 198   switch (c)
 199     {
 200     case 0: return '-';
 201     case 1: return 'H';
 202     case 2: return 'E';
 203     case 3: return 'C';
 204     case 4: return 'S';
 205     case 5: return 'T';
 206     case 6: return 'G';
 207     case 7: return 'B';
 208     case 8: return 'I';
 209     }
 210   return '?';
 211 }
 212
 213
 214 //////////////////////////////////////////////////////////////////////////////
 215 /**
 216  * @brief Transforms the solvend accessiblity code into an integer number
 217  */
 218 inline char
 219 sa2i(char c)
 220 {
 221   //- A B C D E
 222   if (c>='a' && c<='z') c+='A'-'a';
 223   switch (c)
 224     {
 225     case '.': return 0;
 226     case '-': return 0;
 227     case 'A': return 1;
 228     case 'B': return 2;
 229     case 'C': return 3;
 230     case 'D': return 4;
 231     case 'E': return 5;
 232     case 'F': return 6;
 233     case ' ': return -1;
 234     case '\t': return -1;
 235     case '\n': return -1;
 236     }
 237   return -2;
 238 }
 239
 240 //////////////////////////////////////////////////////////////////////////////
 241 /**
 242  * @brief Transforms integers between 0 and 5 into the solvent accessibility code
 243  */
 244 inline char i2sa(int c)
 245 {
 246   //- H E C S T G B
 247   switch (c)
 248     {
 249     case 0: return '-';
 250     case 1: return 'A';
 251     case 2: return 'B';
 252     case 3: return 'C';
 253     case 4: return 'D';
 254     case 5: return 'E';
 255     case 6: return 'F';
 256     }
 257   return '?';
 258 }
 259
 260
 261 //////////////////////////////////////////////////////////////////////////////
 262 /**
 263  * @brief Transforms alternative secondary structure symbols into symbols
 264  */
 265 inline char
 266 ss2ss(char c)
 267 {
 268   //- H E C S T G B
 269   switch (c)
 270     {
 271     case '~': return 'C';
 272     case 'I': return 'C';
 273     case 'i': return 'c';
 274     case 'H':
 275     case 'E':
 276     case 'C':
 277     case 'S':
 278     case 'T':
 279     case 'G':
 280     case 'B':
 281     case 'h':
 282     case 'e':
 283     case 'c':
 284     case 's':
 285     case 't':
 286     case 'g':
 287     case 'b':
 288     case '.':
 289       return c;
 290     }
 291   return '-';
 292 }
 293
 294 //////////////////////////////////////////////////////////////////////////////
 295 /**
 296  * @brief Transforms confidence values of psipred into internal code
 297  */
 298 inline char
 299 cf2i(char c)
 300 {
 301   switch (c)
 302     {
 303     case '-': return 0;
 304     case '.': return 0;
 305     case '0': return 1;
 306     case '1': return 2;
 307     case '2': return 3;
 308     case '3': return 4;
 309     case '4': return 5;
 310     case '5': return 6;
 311     case '6': return 7;
 312     case '7': return 8;
 313     case '8': return 9;
 314     case '9': return 10;
 315     }
 316   return 0;
 317 }
 318
 319 //////////////////////////////////////////////////////////////////////////////
 320 /**
 321  * @brief Transforms internal representation of psipred confidence values into printable chars
 322  */
 323 inline char
 324 i2cf(char c)
 325 {
 326   switch (c)
 327     {
 328     case 0: return '-';
 329     case 1: return '0';
 330     case 2: return '1';
 331     case 3: return '2';
 332     case 4: return '3';
 333     case 5: return '4';
 334     case 6: return '5';
 335     case 7: return '6';
 336     case 8: return '7';
 337     case 9: return '8';
 338     case 10: return '9';
 339     }
 340   return '-';
 341 }
 342
 343
 344 //////////////////////////////////////////////////////////////////////////////
 345 /**
 346  * @brief Fast lookup of log2(1+2^(-x)) for x>=0 (precision < 0.35%)
 347  */
 348 inline float
 349 fast_addscore(float x)
 350 {
 351   static float val[2001];         // val[i]=log2(1+2^(-x))
 352   static char initialized;
 353   if (x>20) return 0.0;
 354   if (x<0)
 355     {
 356       fprintf(stderr,"Error in function fast_addscore: argument %g is negative\n",x);
 357       exit(7);
 358     }
 359   if (!initialized)   //First fill in the log2-vector
 360     {
 361       for (int i=0; i<=2000; i++) val[i]=log2(1.0+pow(2,-0.01*(i+0.5)));
 362       initialized=1;
 363     }
 364   return val[(int)(100.0*x)];
 365 }
 366
 367
 368
 369 //////////////////////////////////////////////////////////////////////////////
 370 /**
 371  * @brief Little utilities for output
 372  */
 373 inline void
 374 fout(FILE* outf, int d)
 375 {
 376   if (d>=99999) fprintf(outf,"*\t"); else fprintf(outf,"%i\t",d);
 377   return;
 378 }
 379
 380 //////////////////////////////////////////////////////////////////////////////
 381 /**
 382  * @brief Errors
 383  */
 384 int
 385 FormatError(const char infile[], const char details[]="")
 386 {
 387   cerr<<"Error in "<</*par.argv[0],FS*/__FILE__<<": wrong format while reading file \'"<<infile<<". "<<details<<"\n";
 388   exit(1);
 389 }
 390
 391 int
 392 OpenFileError(const char outfile[])
 393 {
 394   cerr<<endl<<"Error in "<</*par.argv[0],FS*/__FILE__<<": could not open file \'"<<outfile<<"\'\n";
 395   exit(2);
 396 }
 397
 398 int
 399 MemoryError(const char arrayname[])
 400 {
 401   cerr<<"Error in "<</*par.argv[0],FS*/__FILE__<<": Memory overflow while creating \'"<<arrayname<<"\'. Please report this bug to developers\n";
 402   exit(3);
 403 }
 404
 405 int
 406 InternalError(const char errstr[])
 407 {
 408   cerr<<"Error in "<</*par.argv[0],FS*/__FILE__<<":  "<<errstr<<". Please report this bug to developers\n";
 409   exit(6);
 410 }
 411
 412
 413 //////////////////////////////////////////////////////////////////////////////
 414 /**
 415  * @brief Takes family code (eg. a.1.2.3) and returns strings 'a', 'a.1', and 'a.1.2'
 416  */
 417 inline void
 418 ScopID(char cl[], char fold[], char sfam[], const char fam[])
 419 {
 420   char* ptr;
 421
 422   //get scop class ID
 423   strcpy(cl,fam);
 424   ptr = strchr(cl,'.');               //return adress of next '.' in name
 425   if(ptr) ptr[0]='\0';
 426
 427   //get scop fold ID
 428   strcpy(fold,fam);
 429   ptr = strchr(fold,'.');             //return adress of next '.' in name
 430   if(ptr) ptr = strchr(ptr+1,'.');    //return adress of next '.' in name
 431   if(ptr) ptr[0]='\0';
 432
 433   //get scop superfamily ID
 434   strcpy(sfam,fam);
 435   ptr = strchr(sfam,'.');            //return adress of next '.' in name
 436   if(ptr) ptr = strchr(ptr+1,'.');   //return adress of next '.' in name
 437   if(ptr) ptr = strchr(ptr+1,'.');   //return adress of next '.' in name
 438   if(ptr) ptr[0]='\0';
 439   return;
 440 }
 441
 442 //////////////////////////////////////////////////////////////////////////////
 443 /**
 444  * @brief Read up to n lines of outfile and write to screen (STDERR)
 445  */
 446 void
 447 WriteToScreen(char* outfile, int n)
 448 {
 449   char line[LINELEN]="";
 450   ifstream outf;
 451   outf.open(outfile, ios::in);
 452   if (!outf) {OpenFileError(outfile);}
 453   cout<<"\n";
 454   for(; n>0 && outf.getline(line,LINELEN); n--) cout<<line<<"\n";
 455   outf.close();
 456   cout<<"\n";
 457 }
 458
 459 inline void
 460 WriteToScreen(char* outfile) {WriteToScreen(outfile,INT_MAX);}
 461
 462
 463
 464 /////////////////////////////////////////////////////////////////////////////////////
 465 /**
 466  * @brief Read .hhdefaults file into array argv_conf (beginning at argv_conf[1])
 467  */
 468 void
 469 ReadDefaultsFile(int& argc_conf, char** argv_conf)
 470 {
 471   char line[LINELEN]="";
 472   char filename[NAMELEN];
 473   char* c_first;   //pointer to first character of argument string
 474   char* c;         //pointer to scan line read in for end of argument
 475   //  ifstream configf;
 476   FILE* configf=NULL;
 477   argc_conf=1;     //counts number of arguments read in
 478
 479   // Open config file
 480   strcpy(filename,"./.hhdefaults");
 481   configf = fopen(filename,"r");
 482   if (!configf && getenv("HOME"))
 483     {
 484       strcpy(filename,getenv("HOME"));
 485       strcat(filename,"/.hhdefaults");
 486       configf = fopen(filename,"r");
 487       if (!configf)
 488         {
 489           if (v>=3) cerr<<"Warning: could not find ./.hhdefaults or "<<filename<<"\n";
 490           return;
 491         }
 492     }
 493   else if (!configf) return; // only webserver has no home directory => need no warning
 494
 495   // Scan file until line 'program_nameANYTHING'
 496   while (fgets(line,LINELEN,configf))
 497     if (!strncmp(line,program_name,6)) break;
 498   // Found line 'program_nameANYTHING'?
 499   if (!strncmp(line,program_name,6))
 500     {
 501       // Read in options until end-of-file or empty line
 502         while (fgets(line,LINELEN,configf) && strcmp(line,"\n"))
 503         {
 504           // Analyze line
 505           c=line;
 506           do
 507             {
 508               // Find next word
 509               while (*c==' ' || *c=='\t') c++; //Advance until next non-white space
 510               if (*c=='\0' || *c=='\n' || *c=='#') break;  //Is next word empty string?
 511               c_first=c;
 512               while (*c!=' ' && *c!='\t'  && *c!='#' && *c!='\0' && *c!='\n' ) c++; //Advance until next white space or '#'
 513               if (*c=='\0' || *c=='\n' || *c=='#')         //Is end of line reached?
 514                 {
 515                   *c='\0';
 516                   argv_conf[argc_conf]=new char[strlen(c_first)+1];
 517                   strcpy(argv_conf[argc_conf++],c_first);
 518                   break;
 519                 }
 520               *c='\0';
 521               argv_conf[argc_conf]=new char[strlen(c_first)+1];
 522               strcpy(argv_conf[argc_conf++],c_first);
 523               printf("Argument: %s\n",c_first);
 524               c++;
 525             } while (1);
 526         } //end read line
 527      if (v>=3)
 528         {
 529           cout<<"Arguments read in from .hhdefaults:";
 530           for (int argc=1; argc<argc_conf; argc++) cout<<(argv_conf[argc][0]=='-'? " ":"")<<argv_conf[argc]<<" ";
 531           cout<<"\n";
 532         }
 533      else if (v>=3) cout<<"Read in "<<argc_conf<<" default arguments for "<<program_name<<" from "<<filename<<"\n";
 534      }
 535   else //found no line 'program_name   anything"
 536     {
 537       if (v>=3) cerr<<endl<<"Warning: no default options for \'"<<program_name<<"\' found in "<<filename<<"\n";
 538       return; //no line 'program_name   anything' found
 539     }
 540 //   configf.close();
 541   fclose(configf);
 542 }
 543
 544
 545 /////////////////////////////////////////////////////////////////////////////////////
 546 /**
 547  * @brief Set default parameter values
 548  */
 549 void
 550 SetDefaults(hhalign_para rHhalignPara)
 551 {
 552
 553   par.append=0;                // overwrite output file
 554   par.outformat=0;             // 0: hhr  1: FASTA  2:A2M   3:A3M
 555   par.p=20.0f;                 // minimum threshold for inclusion in hit list and alignment listing
 556   par.E=1e6f;                  // maximum threshold for inclusion in hit list and alignment listing
 557   par.b=10;                    // min number of alignments
 558   par.B=500;                   // max number of alignments
 559   par.z=10;                    // min number of lines in hit list
 560   par.Z=500;                   // max number of lines in hit list
 561   par.e=1e-3f;                 // maximum E-value for inclusion in output alignment, output HMM, and PSI-BLAST checkpoint model
 562   par.showcons=1;              // show consensus sequence
 563   par.showdssp=1;              // show predicted secondary structure
 564   par.showpred=1;              // show dssp secondary structure
 565   par.cons=0;                  // show first non-SS sequence as main representative sequence (not consensus)
 566   par.nseqdis=1;               // maximum number of query sequences for output alignment
 567   par.mark=0;                  // 1: only marked sequences (or first) get displayed; 0: most divergent ones get displayed
 568   par.aliwidth=80;             // number of characters per line in output alignments for HMM search
 569   par.max_seqid=90;            // default for maximum sequence identity threshold
 570   par.qid=0;                   // default for minimum sequence identity with query
 571   par.qsc=-20.0f;              // default for minimum score per column with query
 572   par.coverage=0;              // default for minimum coverage threshold
 573   par.Ndiff=100;               // pick Ndiff most different sequences from alignment
 574   par.coverage_core=80;        // Minimum coverage for sequences in core alignment
 575   par.qsc_core=0.3f;           // Minimum score per column of core sequence with query
 576   par.coresc=-20.0f;           // Minimum score per column with core alignment (HMM)
 577
 578   par.M=1;                     // match state assignment is by A2M/A3M
 579   par.Mgaps=50;                // Above this percentage of gaps, columns are assigned to insert states (for par.M=2)
 580   par.calibrate=0;             // default: no calibration
 581   par.calm=0;                  // derive P-values from: 0:query calibration  1:template calibration  2:both
 582   par.mode=0;                  //
 583
 584   par.wg=0;                    // 0: use local sequence weights   1: use local ones
 585
 586   par.matrix=0;                // Subst.matrix 0: Gonnet, 1: HSDM, 2: BLOSUM50 3: BLOSUM62
 587   par.pcm=2;                   // pseudocount mode: default=divergence-dependent (but not column-specific)
 588 #if 1 /* Nelder-Meade on Baliscore*/
 589   par.pca=1.712190f;                // default values for substitution matrix pseudocounts
 590   par.pcb=1.039640f;                // significant reduction of pcs by Neff_M starts around Neff_M-1=pcb
 591   par.pcc=0.878067f;                // pcs are reduced prop. to 1/Neff^pcc
 592   par.pcw=0.0f;                // wc>0 weighs columns according to their intra-clomun similarity
 593
 594   par.gapb=1.405220;                // default values for transition pseudocounts
 595   par.gapd=1.316760;               // gap open penalty pseudocount; 0.25 corresponds to 7.1*gapf bits
 596   par.gape=1.793780;                // gap extension penalty pseudocount
 597   par.gapf=1.034710;                // factor for increasing gap open penalty for deletes
 598   par.gapg=0.894772;                // factor for increasing gap open penalty for inserts
 599   par.gaph=0.544072;                // factor for increasing gap extension penalty for deletes
 600   par.gapi=0.862559;                // factor for increasing gap extension penalty for inserts
 601 #else  /* Soeding's default*/
 602   par.pca=1.0f;                // default values for substitution matrix pseudocounts
 603   par.pcb=1.5f;                // significant reduction of pcs by Neff_M starts around Neff_M-1=pcb
 604   par.pcc=1.0f;                // pcs are reduced prop. to 1/Neff^pcc
 605   par.pcw=0.0f;                // wc>0 weighs columns according to their intra-clomun similarity
 606
 607   par.gapb=1.0;                // default values for transition pseudocounts
 608   par.gapd=0.15;               // gap open penalty pseudocount; 0.25 corresponds to 7.1*gapf bits
 609   par.gape=1.0;                // gap extension penalty pseudocount
 610   par.gapf=0.6;                // factor for increasing gap open penalty for deletes
 611   par.gapg=0.6;                // factor for increasing gap open penalty for inserts
 612   par.gaph=0.6;                // factor for increasing gap extension penalty for deletes
 613   par.gapi=0.6;                // factor for increasing gap extension penalty for inserts
 614 #endif
 615
 616 #if 0
 617   /* Viterbi parameters optimised on Sabre (R228), FS, r228 -> r229 */
 618   par.pcaV=1.245150f;                // default values for substitution matrix pseudocounts
 619   par.pcbV=1.682110f;                // significant reduction of pcs by Neff_M starts around Neff_M-1=pcb
 620   par.pccV=1.483840f;                // pcs are reduced prop. to 1/Neff^pcc
 621   par.pcwV=0.0f;                // wc>0 weighs columns according to their intra-clomun similarity
 622
 623   par.gapbV=0.818625;                // default values for transition pseudocounts
 624   par.gapdV=0.666110;               // gap open penalty pseudocount; 0.25 corresponds to 7.1*gapf bits
 625   par.gapeV=1.028050;                // gap extension penalty pseudocount
 626   par.gapfV=0.710760;                // factor for increasing gap open penalty for deletes
 627   par.gapgV=1.649800;                // factor for increasing gap open penalty for inserts
 628   par.gaphV=0.470604;                // factor for increasing gap extension penalty for deletes
 629   par.gapiV=0.829479;                // factor for increasing gap extension penalty for inserts
 630 #elif 1
 631   /* Viterbi parameters optimised on Balibase, r244 -> r245 */
 632   par.pcaV=1.333860f;                // default values for substitution matrix pseudocounts
 633   par.pcbV=1.934480f;                // significant reduction of pcs by Neff_M starts around Neff_M-1=pcb
 634   par.pccV=1.655610f;                // pcs are reduced prop. to 1/Neff^pcc
 635   par.pcwV=0.0f;                // wc>0 weighs columns according to their intra-clomun similarity
 636
 637   par.gapbV=0.334525;                // default values for transition pseudocounts
 638   par.gapdV=0.074534;               // gap open penalty pseudocount; 0.25 corresponds to 7.1*gapf bits
 639   par.gapeV=0.320336;                // gap extension penalty pseudocount
 640   par.gapfV=0.151634;                // factor for increasing gap open penalty for deletes
 641   par.gapgV=0.641516;                // factor for increasing gap open penalty for inserts
 642   par.gaphV=0.266434;                // factor for increasing gap extension penalty for deletes
 643   par.gapiV=0.598414;                // factor for increasing gap extension penalty for inserts
 644 #else /* Soeding default*/
 645   par.pcaV=1.0f;                // default values for substitution matrix pseudocounts
 646   par.pcbV=1.5f;                // significant reduction of pcs by Neff_M starts around Neff_M-1=pcb
 647   par.pccV=1.0f;                // pcs are reduced prop. to 1/Neff^pcc
 648   par.pcwV=0.0f;                // wc>0 weighs columns according to their intra-clomun similarity
 649
 650   par.gapbV=1.0;                // default values for transition pseudocounts
 651   par.gapdV=0.15;               // gap open penalty pseudocount; 0.25 corresponds to 7.1*gapf bits
 652   par.gapeV=1.0;                // gap extension penalty pseudocount
 653   par.gapfV=0.6;                // factor for increasing gap open penalty for deletes
 654   par.gapgV=0.6;                // factor for increasing gap open penalty for inserts
 655   par.gaphV=0.6;                // factor for increasing gap extension penalty for deletes
 656   par.gapiV=0.6;                // factor for increasing gap extension penalty for inserts
 657 #endif
 658
 659   par.ssm=2;                   // ss scoring mode: 0:no ss score  1:score after alignment  2:score during alignment
 660   par.ssw=0.11f;               // weight of ss scoring
 661   par.ssa=1.0f;                // weight of ss evolution matrix
 662   par.shift=-0.01f;            // Shift match score up
 663   par.mact=0.3001f;            // Score threshold for MAC alignment in local mode (set to 0.5001 to track user modification)
 664   par.corr=0.1f;               // Weight of correlations of scores for |i-j|<=4
 665   par.wstruc=1.0f;             // Weight of structure scores
 666
 667   par.egq=0.0f;                // no charge for end gaps as default
 668   par.egt=0.0f;                // no charge for end gaps as default
 669
 670   par.trans=0;                 // no transitive scoring as default
 671   par.Emax_trans=100.0f;       // use intermediate HMMs with E-values up to 100 between query and database HMM
 672   par.Emax_trans=100.0f;       // use intermediate HMMs with E-values up to 100 between query and database HMM
 673   par.wtrans=1.0f;             // Ztot[k] = Zq[k] + wtrans * (Zforward[k]+Zreverse[k])
 674   par.ssgap=0;                 // 1: add secondary structure-dependent gap penalties  0:off
 675   par.ssgapd=1.0f;             // secondary structure-dependent gap-opening penalty (per residue)
 676   par.ssgape=0.0f;             // secondary structure-dependent gap-extension penalty (per residue)
 677   par.ssgapi=4;                // max. number of inside-integer(ii); gap-open-penalty= -ii*ssgapd
 678
 679   par.loc=1;                   // local vs. global alignment as default
 680   par.altali=2;                // find up to two (possibly overlapping) subalignments
 681   par.forward=0;               // 0: Viterbi algorithm; 1: Viterbi+stochastic sampling; 3:Maximum Accuracy (MAC) algorithm
 682   par.realign=1;               // realign with MAC algorithm
 683
 684   par.repmode=0;               // repeats score independently of one another
 685   par.columnscore=1;           // Default column score is 1: null model pnul = 1/2 * (q_av(a)+p_av(a))
 686   par.min_overlap=0;           // automatic minimum overlap used
 687   par.opt=0;                   // Default = optimization mode off
 688   par.readdefaultsfile=0;      // Default = do not read a defaults file ./.hhdefaults or HOME/.hhdefaults
 689   par.maxdbstrlen=200;         // maximum length of database string to be printed in 'Command' line of hhr file
 690   par.mode=0;
 691   par.idummy=par.jdummy=0;     //
 692
 693   par.notags=1;                // neutralize His-tags, FLAG-tags, C-myc-tags
 694
 695   // Initialize strings
 696   strcpy(par.infile,"stdin");
 697   strcpy(par.outfile,"");
 698   strcpy(par. pairwisealisfile,"");
 699   strcpy(par.buffer,"buffer.txt");
 700   strcpy(par.scorefile,"");
 701   strcpy(par.wfile,"");
 702   strcpy(par.alnfile,"");
 703   strcpy(par.hhmfile,"");
 704   strcpy(par.psifile,"");
 705   par.exclstr=NULL;
 706
 707
 708   //#if PARAMETERSFROMFILE  /* read parameter file from home-dir */
 709   //#include "hhutil-C-help.h"
 710   //#endif /* read parameter file from home-dir */
 711   if (rHhalignPara.pca  >= 0.00){ par.pca  = rHhalignPara.pca;  }
 712   if (rHhalignPara.pcb  >= 0.00){ par.pcb  = rHhalignPara.pcb;  }
 713   if (rHhalignPara.pcc  >= 0.00){ par.pcc  = rHhalignPara.pcc;  }
 714   if (rHhalignPara.pcw  >= 0.00){ par.pcw  = rHhalignPara.pcw;  }
 715   if (rHhalignPara.gapb >= 0.00){ par.gapb = rHhalignPara.gapb; }
 716   if (rHhalignPara.gapd >= 0.00){ par.gapd = rHhalignPara.gapd; }
 717   if (rHhalignPara.gape >= 0.00){ par.gape = rHhalignPara.gape; }
 718   if (rHhalignPara.gapf >= 0.00){ par.gapf = rHhalignPara.gapf; }
 719   if (rHhalignPara.gapg >= 0.00){ par.gapg = rHhalignPara.gapg; }
 720   if (rHhalignPara.gaph >= 0.00){ par.gaph = rHhalignPara.gaph; }
 721   if (rHhalignPara.gapi >= 0.00){ par.gapi = rHhalignPara.gapi; }
 722   if (rHhalignPara.pcaV  >= 0.00){ par.pcaV  = rHhalignPara.pcaV;  }
 723   if (rHhalignPara.pcbV  >= 0.00){ par.pcbV  = rHhalignPara.pcbV;  }
 724   if (rHhalignPara.pccV  >= 0.00){ par.pccV  = rHhalignPara.pccV;  }
 725   if (rHhalignPara.pcwV  >= 0.00){ par.pcwV  = rHhalignPara.pcwV;  }
 726   if (rHhalignPara.gapbV >= 0.00){ par.gapbV = rHhalignPara.gapbV; }
 727   if (rHhalignPara.gapdV >= 0.00){ par.gapdV = rHhalignPara.gapdV; }
 728   if (rHhalignPara.gapeV >= 0.00){ par.gapeV = rHhalignPara.gapeV; }
 729   if (rHhalignPara.gapfV >= 0.00){ par.gapfV = rHhalignPara.gapfV; }
 730   if (rHhalignPara.gapgV >= 0.00){ par.gapgV = rHhalignPara.gapgV; }
 731   if (rHhalignPara.gaphV >= 0.00){ par.gaphV = rHhalignPara.gaphV; }
 732   if (rHhalignPara.gapiV >= 0.00){ par.gapiV = rHhalignPara.gapiV; }
 733
 734    return;
 735 } /** this is the end of SetDefaults() **/
 736
 737 void SetRnaDefaults(hhalign_para rHhalignPara)
 738 {
 739   par.pca=1.28f;                // default values for substitution matrix pseudocounts
 740   par.pcb=1.75f;                // significant reduction of pcs by Neff_M starts around Neff_M-1=pcb
 741   par.pcc=0.87f;                // pcs are reduced prop. to 1/Neff^pcc
 742   par.pcw=0.0f;                // wc>0 weighs columns according to their intra-clomun similarity
 743
 744   par.gapb=0.80;                // default values for transition pseudocounts
 745   par.gapd=0.34;               // gap open penalty pseudocount; 0.25 corresponds to 7.1*gapf bits
 746   par.gape=2.25;                // gap extension penalty pseudocount
 747   par.gapf=0.51;                // factor for increasing gap open penalty for deletes
 748   par.gapg=1.01;                // factor for increasing gap open penalty for inserts
 749   par.gaph=1.24;                // factor for increasing gap extension penalty for deletes
 750   par.gapi=0.45;                // factor for increasing gap extension penalty for inserts
 751
 752 #if 0  /* these are the parameters determined by Dave (pre r274) */
 753   par.pcaV=2.57f;                // default values for substitution matrix pseudocounts
 754   par.pcbV=2.34f;                // significant reduction of pcs by Neff_M starts around Neff_M-1=pcb
 755   par.pccV=0.88f;                // pcs are reduced prop. to 1/Neff^pcc
 756   par.pcwV=0.0f;                // wc>0 weighs columns according to their intra-clomun similarity
 757
 758   par.gapbV=1.41;                // default values for transition pseudocounts
 759   par.gapdV=1.8;               // gap open penalty pseudocount; 0.25 corresponds to 7.1*gapf bits
 760   par.gapeV=1.5;                // gap extension penalty pseudocount
 761   par.gapfV=1.03;                // factor for increasing gap open penalty for deletes
 762   par.gapgV=0.89;                // factor for increasing gap open penalty for inserts
 763   par.gaphV=0.54;                // factor for increasing gap extension penalty for deletes
 764   par.gapiV=0.86;                // factor for increasing gap extension penalty for inserts
 765 #else /* parameters determined for r274, using Bralibase, FS */
 766   par.pcaV=1.655620f;                // default values for substitution matrix pseudocounts
 767   par.pcbV=0.438399f;                // significant reduction of pcs by Neff_M starts around Neff_M-1=pcb
 768   par.pccV=0.371491f;                // pcs are reduced prop. to 1/Neff^pcc
 769   par.pcwV=0.0f;                // wc>0 weighs columns according to their intra-clomun similarity
 770
 771   par.gapbV=1.914490;                // default values for transition pseudocounts
 772   par.gapdV=0.104278;               // gap open penalty pseudocount; 0.25 corresponds to 7.1*gapf bits
 773   par.gapeV=1.100210;                // gap extension penalty pseudocount
 774   par.gapfV=0.335152;                // factor for increasing gap open penalty for deletes
 775   par.gapgV=0.786688;                // factor for increasing gap open penalty for inserts
 776   par.gaphV=0.667143;                // factor for increasing gap extension penalty for deletes
 777   par.gapiV=0.711993;                // factor for increasing gap extension penalty for inserts
 778 #endif
 779
 780   //#if PARAMETERSFROMFILE  /* read parameter file from home-dir */
 781   //#include "hhutil-C-help.h"
 782   //#endif /* read parameter file from home-dir */
 783   if (rHhalignPara.pca  >= 0.00){ par.pca  = rHhalignPara.pca;  }
 784   if (rHhalignPara.pcb  >= 0.00){ par.pcb  = rHhalignPara.pcb;  }
 785   if (rHhalignPara.pcc  >= 0.00){ par.pcc  = rHhalignPara.pcc;  }
 786   if (rHhalignPara.pcw  >= 0.00){ par.pcw  = rHhalignPara.pcw;  }
 787   if (rHhalignPara.gapb >= 0.00){ par.gapb = rHhalignPara.gapb; }
 788   if (rHhalignPara.gapd >= 0.00){ par.gapd = rHhalignPara.gapd; }
 789   if (rHhalignPara.gape >= 0.00){ par.gape = rHhalignPara.gape; }
 790   if (rHhalignPara.gapf >= 0.00){ par.gapf = rHhalignPara.gapf; }
 791   if (rHhalignPara.gapg >= 0.00){ par.gapg = rHhalignPara.gapg; }
 792   if (rHhalignPara.gaph >= 0.00){ par.gaph = rHhalignPara.gaph; }
 793   if (rHhalignPara.gapi >= 0.00){ par.gapi = rHhalignPara.gapi; }
 794   if (rHhalignPara.pcaV  >= 0.00){ par.pcaV  = rHhalignPara.pcaV;  }
 795   if (rHhalignPara.pcbV  >= 0.00){ par.pcbV  = rHhalignPara.pcbV;  }
 796   if (rHhalignPara.pccV  >= 0.00){ par.pccV  = rHhalignPara.pccV;  }
 797   if (rHhalignPara.pcwV  >= 0.00){ par.pcwV  = rHhalignPara.pcwV;  }
 798   if (rHhalignPara.gapbV >= 0.00){ par.gapbV = rHhalignPara.gapbV; }
 799   if (rHhalignPara.gapdV >= 0.00){ par.gapdV = rHhalignPara.gapdV; }
 800   if (rHhalignPara.gapeV >= 0.00){ par.gapeV = rHhalignPara.gapeV; }
 801   if (rHhalignPara.gapfV >= 0.00){ par.gapfV = rHhalignPara.gapfV; }
 802   if (rHhalignPara.gapgV >= 0.00){ par.gapgV = rHhalignPara.gapgV; }
 803   if (rHhalignPara.gaphV >= 0.00){ par.gaphV = rHhalignPara.gaphV; }
 804   if (rHhalignPara.gapiV >= 0.00){ par.gapiV = rHhalignPara.gapiV; }
 805
 806 } /* this is the end of SetRnaDefaults() */
 807
 808 void SetDnaDefaults(hhalign_para rHhalignPara)
 809 {
 810   par.pca=2.89f;                // default values for substitution matrix pseudocounts
 811   par.pcb=1.17f;                // significant reduction of pcs by Neff_M starts around Neff_M-1=pcb
 812   par.pcc=0.88f;                // pcs are reduced prop. to 1/Neff^pcc
 813   par.pcw=0.0f;                // wc>0 weighs columns according to their intra-clomun similarity
 814
 815   par.gapb=0.80;                // default values for transition pseudocounts
 816   par.gapd=0.34;               // gap open penalty pseudocount; 0.25 corresponds to 7.1*gapf bits
 817   par.gape=2.25;                // gap extension penalty pseudocount
 818   par.gapf=0.78;                // factor for increasing gap open penalty for deletes
 819   par.gapg=1.01;                // factor for increasing gap open penalty for inserts
 820   par.gaph=1.24;                // factor for increasing gap extension penalty for deletes
 821   par.gapi=0.45;                // factor for increasing gap extension penalty for inserts
 822
 823 #if 0 /* these are the parameters determined by Dave (pre r274) */
 824   par.pcaV=1.712f;                // default values for substitution matrix pseudocounts
 825   par.pcbV=1.039f;                // significant reduction of pcs by Neff_M starts around Neff_M-1=pcb
 826   par.pccV=0.266f;                // pcs are reduced prop. to 1/Neff^pcc
 827   par.pcwV=0.0f;                // wc>0 weighs columns according to their intra-clomun similarity
 828
 829   par.gapbV=1.405;                // default values for transition pseudocounts
 830   par.gapdV=1.8;               // gap open penalty pseudocount; 0.25 corresponds to 7.1*gapf bits
 831   par.gapeV=2.25;                // gap extension penalty pseudocount
 832   par.gapfV=1.034;                // factor for increasing gap open penalty for deletes
 833   par.gapgV=2.025;                // factor for increasing gap open penalty for inserts
 834   par.gaphV=0.544;                // factor for increasing gap extension penalty for deletes
 835   par.gapiV=1.35;                // factor for increasing gap extension penalty for inserts }
 836 #else /* parameters determined for r274, using mdsa, FS */
 837   par.pcaV=2.196;                // default values for substitution matrix pseudocounts
 838   par.pcbV=0.329;                // significant reduction of pcs by Neff_M starts around Neff_M-1=pcb
 839   par.pccV=0.393;                // pcs are reduced prop. to 1/Neff^pcc
 840   par.pcwV=0.0f;                // wc>0 weighs columns according to their intra-clomun similarity
 841
 842   par.gapbV=0.570;                // default values for transition pseudocounts
 843   par.gapdV=0.048;               // gap open penalty pseudocount; 0.25 corresponds to 7.1*gapf bits
 844   par.gapeV=1.692;                // gap extension penalty pseudocount
 845   par.gapfV=0.398;                // factor for increasing gap open penalty for deletes
 846   par.gapgV=0.121;                // factor for increasing gap open penalty for inserts
 847   par.gaphV=0.012;                // factor for increasing gap extension penalty for deletes
 848   par.gapiV=0.645;                // factor for increasing gap extension penalty for inserts }
 849 #endif
 850
 851   //#if PARAMETERSFROMFILE  /* read parameter file from home-dir */
 852   //#include "hhutil-C-help.h"
 853   //#endif /* read parameter file from home-dir */
 854   if (rHhalignPara.pca  >= 0.00){ par.pca  = rHhalignPara.pca;  }
 855   if (rHhalignPara.pcb  >= 0.00){ par.pcb  = rHhalignPara.pcb;  }
 856   if (rHhalignPara.pcc  >= 0.00){ par.pcc  = rHhalignPara.pcc;  }
 857   if (rHhalignPara.pcw  >= 0.00){ par.pcw  = rHhalignPara.pcw;  }
 858   if (rHhalignPara.gapb >= 0.00){ par.gapb = rHhalignPara.gapb; }
 859   if (rHhalignPara.gapd >= 0.00){ par.gapd = rHhalignPara.gapd; }
 860   if (rHhalignPara.gape >= 0.00){ par.gape = rHhalignPara.gape; }
 861   if (rHhalignPara.gapf >= 0.00){ par.gapf = rHhalignPara.gapf; }
 862   if (rHhalignPara.gapg >= 0.00){ par.gapg = rHhalignPara.gapg; }
 863   if (rHhalignPara.gaph >= 0.00){ par.gaph = rHhalignPara.gaph; }
 864   if (rHhalignPara.gapi >= 0.00){ par.gapi = rHhalignPara.gapi; }
 865   if (rHhalignPara.pcaV  >= 0.00){ par.pcaV  = rHhalignPara.pcaV;  }
 866   if (rHhalignPara.pcbV  >= 0.00){ par.pcbV  = rHhalignPara.pcbV;  }
 867   if (rHhalignPara.pccV  >= 0.00){ par.pccV  = rHhalignPara.pccV;  }
 868   if (rHhalignPara.pcwV  >= 0.00){ par.pcwV  = rHhalignPara.pcwV;  }
 869   if (rHhalignPara.gapbV >= 0.00){ par.gapbV = rHhalignPara.gapbV; }
 870   if (rHhalignPara.gapdV >= 0.00){ par.gapdV = rHhalignPara.gapdV; }
 871   if (rHhalignPara.gapeV >= 0.00){ par.gapeV = rHhalignPara.gapeV; }
 872   if (rHhalignPara.gapfV >= 0.00){ par.gapfV = rHhalignPara.gapfV; }
 873   if (rHhalignPara.gapgV >= 0.00){ par.gapgV = rHhalignPara.gapgV; }
 874   if (rHhalignPara.gaphV >= 0.00){ par.gaphV = rHhalignPara.gaphV; }
 875   if (rHhalignPara.gapiV >= 0.00){ par.gapiV = rHhalignPara.gapiV; }
 876
 877
 878 } /* this is the end of SetDnaDefaults() */
 879
 880
 881 /*
 882  * EOF hhutil-C.h
 883  */