sources/jnet/jnet.c~

   1 /************************************************************************
   2  *               Jnet - A consensus neural network protein              *
   3  *                      secondary structure prediction method           *
   4  *                                                                      *
   5  *  Copyright 1999,2009 James Cuff, Jonathan Barber and Christian Cole  *
   6  *                                                                      *
   7  ************************************************************************
   8
   9
  10 -------------------------------------------------------------------------
  11     This program is free software: you can redistribute it and/or modify
  12     it under the terms of the GNU General Public License as published by
  13     the Free Software Foundation, either version 3 of the License, or
  14     (at your option) any later version.
  15
  16     This program is distributed in the hope that it will be useful,
  17     but WITHOUT ANY WARRANTY; without even the implied warranty of
  18     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  19     GNU General Public License for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with this program.  If not, see <http://www.gnu.org/licenses/>.
  23 -------------------------------------------------------------------------*/
  24
  25 /* Fri Nov 15 15:28:09 GMT 2002 Bug fixes by JDB */
  26 /* WARNING: MUST BE PRESENTED WITH UPPERCASE DATA as lowercase 'n'
  27  * has a special internal meaning in seq2int
  28  * */
  29
  30                                                                                            /*#define POST *//* Allow post processing of NN output */
  31                                                                                                    /*#define FILTER *//* Allow string post processing of final output */
  32
  33 /** The profile-specific defines below are not independent of each other.
  34 Especially check that the Jury set-up is correct! */
  35 #define HMM                     /* Do HMMer predictions */
  36 #define PSSM                    /* Do PSSM predicitons */
  37                                                                                                    /*#define HMMONLY *//* For non-jury positions use the HMM prediction - only works when DOJURY is undefined */
  38
  39 /* The allowed value of the integers representing residues */
  40 enum
  41 {
  42   WINAR = 30,                   /* array size for winar */
  43   PSILEN = 20,                  /* no. of elements per position in PSIBLAST PSSM */
  44   PROLEN = 24,                  /* no. of elements per position for other profiles */
  45   HELIX = 1,                    /* enums for sec. struct. states */
  46   SHEET = 0,
  47   COIL = 2,
  48   MAXSEQNUM = 1000,             /* Maximum number of sequences per alignment */
  49   MAXSEQLEN = 5000,             /* Maximum Sequence Length */
  50   MAXBLOCK = 200                /* Maximum block size for human readable output */
  51 };
  52
  53 #include <stdio.h>
  54 #include <string.h>
  55 #include <stdlib.h>
  56 #include <math.h>
  57 #include <ctype.h>
  58 #include <assert.h>
  59
  60 #include "cmdline/cmdline.h"
  61
  62 #include "pssma1.h"             /* Uses information from data.psimat, (PSSM) - 100 hidden nodes */
  63 #include "pssma2.h"             /* Second network for pssma1 - 100 hidden nodes */
  64 #include "pssmb1.h"             /* Uses information from data.psimat, (PSSM)  - 20 hidden nodes */
  65 #include "pssmb2.h"             /* Second network for pssmb1  - 20 hidden nodes */
  66
  67 #include "hmm1.h"               /* First network for HMMER profile prediction */
  68 #include "hmm2.h"               /* Second network for hmm1 */
  69
  70 #include "hmmsol25.h"
  71 #include "hmmsol5.h"
  72 #include "hmmsol0.h"
  73
  74 #include "psisol25.h"
  75 #include "psisol5.h"
  76 #include "psisol0.h"
  77
  78 /* Main data structure for information about input data. Not all the fields are
  79  * used */
  80 typedef struct alldata
  81 {
  82   int lens;                     /* length of sequences */
  83   int profmat[MAXSEQLEN][PROLEN];       /* BLOSUM frequencies from PSIBLAST alignment - unused but can't remove without segfault? */
  84   float psimat[MAXSEQLEN][PSILEN];      /* PSIBLAST PSSM */
  85   float hmmmat[MAXSEQLEN][PROLEN];      /* results from HMMER */
  86 }
  87 alldata;
  88
  89 /* Global Vars - used to specify which input profiles have been provided */
  90 int NOPSI = 0;
  91 int NOHMM = 0;
  92 int TESTOUTPUT = 0;
  93 int HUMAN = 0;                  /* toggle human readable output */
  94 int BLOCK = 60;                 /* block size for human output */
  95
  96 /* PSSM-specific functions */
  97 void readpsi (FILE * psifile, alldata * data);
  98 void doprofpsi (const alldata * data, const int win, const int curpos, float psiar[30][30]);
  99
 100 /* HMMer-specific functions */
 101 void readhmm (FILE * hmmfile, alldata * data);
 102 void doprofhmm (const alldata * data, const int win, const int curpos, float psiar[30][30]);
 103
 104 /* General or non-specific functions */
 105 void int2pred (const int *pred, char *letpred, const int length);
 106 float doconf (const float confa, const float confb, const float confc);
 107 void dowinsspsi (float seq2str[MAXSEQLEN][3], int win, int curpos, float winarss[30][4], const alldata * data);
 108 void dopred (const alldata * data);
 109 void test_size (int size, int required_size, char *name);
 110 void print_human (int length, char pred[MAXSEQLEN], float conf[MAXSEQLEN]);
 111
 112 /* Prediction post-processing functions - not be used anymore, but retained for testing */
 113 /* Executing controlled by #define FILTER */
 114 void filter (char inp[MAXSEQLEN], char outp[MAXSEQLEN]);
 115 void StrReplStr (char *s1, char *s2, char *FromSStr, char *ToSStr);
 116 char *Strncat (char *s1, char *s2, int len);
 117
 118 /*************************************************************/
 119 int
 120 main (int argc, char **argv)
 121 {
 122   FILE *psifile, *hmmfile;
 123   alldata *data;
 124
 125   /* Declared in cmdline/cmdline.h */
 126   struct gengetopt_args_info args_info;
 127
 128   /* Set up the data structure for the input data */
 129   if ((data = malloc (sizeof (alldata))) == NULL) {
 130     fprintf (stderr, "ERROR! malloc failed at line %i\n", __LINE__);
 131     exit (EXIT_FAILURE);
 132   }
 133
 134   /* Get the command line arguments */
 135   if (cmdline_parser (argc, argv, &args_info) != 0)
 136     exit (EXIT_FAILURE);
 137
 138   /* check whether TEST STYLE OUTPUT has been requested */
 139   if (args_info.test_flag) {
 140     TESTOUTPUT = 1;
 141   }
 142
 143   /* Print version info - CMDLINE_PARSER_PACKAGE CMDLINE_PARSER_VERSION declared in cmdline/cmdline.h */
 144   if (1 == TESTOUTPUT) {
 145     fprintf (stderr, "%s %s\n", CMDLINE_PARSER_PACKAGE, CMDLINE_PARSER_VERSION);
 146   } else {
 147     fprintf (stdout, "%s %s\n", CMDLINE_PARSER_PACKAGE, CMDLINE_PARSER_VERSION);
 148   }
 149
 150   /* Read the HMMER file. Don't need to check if it is given as is a required option */
 151   if ((hmmfile = fopen (args_info.hmmer_arg, "r")) == NULL) {
 152     fprintf (stderr, "ERROR: Can't open HMMER file %s\n", args_info.hmmer_arg);
 153     exit (EXIT_FAILURE);
 154   } else {
 155     if (1 == TESTOUTPUT) {
 156       fprintf (stderr, "Found HMM profile data\n");
 157     } else {
 158       fprintf (stdout, "Found HMM profile data\n");
 159     }
 160     readhmm (hmmfile, data);
 161     (void) fclose (hmmfile);
 162   }
 163
 164   /* Read the PSIBLAST PSSM */
 165   if (args_info.pssm_given) {
 166     if ((psifile = fopen (args_info.pssm_arg, "r")) == NULL) {
 167       fprintf (stderr, "ERROR: Can't open PSSM profile file %s\n", args_info.pssm_arg);
 168       exit (EXIT_FAILURE);
 169     } else {
 170       if (1 == TESTOUTPUT) {
 171         fprintf (stderr, "Found PSSM profile file\n");
 172       } else {
 173         fprintf (stdout, "Found PSSM profile file\n");
 174       }
 175       readpsi (psifile, data);
 176       (void) fclose (psifile);
 177     }
 178   } else {
 179     NOPSI = 1;
 180   }
 181
 182   /* check whether human readable output is required (set by --human) */
 183   if (args_info.human_given) {
 184     HUMAN = 1;
 185     if (args_info.human_orig != NULL) {
 186       printf ("human: %d\n", args_info.human_arg);
 187       printf ("human orig: %s\n", args_info.human_orig);
 188       BLOCK = args_info.human_arg;
 189       if (BLOCK > MAXBLOCK) {
 190         fprintf (stderr, "ERROR - argument for --human is too large\n");
 191         exit (EXIT_FAILURE);
 192       }
 193     }
 194   }
 195
 196   /* Do the prediction */
 197   if (1 == TESTOUTPUT) {
 198     fprintf (stderr, "Running final predictions!\n");
 199   } else {
 200     fprintf (stdout, "Running final predictions!\n");
 201   }
 202   dopred (data);
 203
 204   free (data);
 205
 206   return EXIT_SUCCESS;
 207 }
 208
 209
 210 /*************************************************************************/
 211 /* Replaces particular runs of secondary structure with something else */
 212 void
 213 filter (char inp[MAXSEQLEN], char outp[MAXSEQLEN])
 214 {
 215   int i;
 216
 217   struct filter_pair
 218   {
 219     char *orig;
 220     char *final;
 221   } pairs[] = {
 222     {
 223     "EHHHE", "EEEEE"}, {
 224     "-HHH-", "HHHHH"}, {
 225     "EHHH-", "EHHHH"}, {
 226     "HHHE-", "HHH--"}, {
 227     "-EHHH", "--HHH"}, {
 228     "EHHE", "EEEE"}, {
 229     "-HEH-", "-HHH-"}, {
 230     "EHH-", "EEE-"}, {
 231     "-HHE", "-EEE"}, {
 232     "-HH-", "----"}, {
 233     "HEEH", "EEEE"}, {
 234     "-HE", "--E"}, {
 235     "EH-", "E--"}, {
 236     "-H-", "---"}, {
 237     "HEH", "HHH"}, {
 238     "-E-E-", "-EEE-"}, {
 239     "E-E", "EEE"}, {
 240     "H-H", "HHH"}, {
 241     "EHE", "EEE"}, {
 242     NULL, NULL}
 243   };
 244
 245   for (i = 0; pairs[i].orig != NULL; i++) {
 246     StrReplStr (outp, inp, pairs[i].orig, pairs[i].final);
 247     strcpy (inp, outp);
 248   }
 249 }
 250
 251 /*************************************************************************/
 252 /* Reads in the PSIBLAST PSSM */
 253 void
 254 readpsi (FILE * psifile, alldata * data)
 255 {
 256   int nlines = 0;
 257
 258   while ((getc (psifile)) != EOF) {
 259     int i;
 260     for (i = 0; i < 20; i++)
 261       fscanf (psifile, "%f", &data->psimat[nlines][i]);
 262     nlines++;
 263   }
 264
 265   /* check the data length matches that of the HMMer data */
 266   if (data->lens != nlines - 1) {
 267     fprintf (stderr, "ERROR! PSSM data length doesn't match that for HMMer. Are both input files from the same alignment/sequence?\n");
 268     exit (EXIT_FAILURE);
 269   }
 270 }
 271
 272
 273 /*************************************************************************/
 274 /* Reads in HMMER profile that was originally based upon ClustalW alignments */
 275 /* Additionally, use this data to define the length of the sequence in the absence of an alignment */
 276 void
 277 readhmm (FILE * hmmfile, alldata * data)
 278 {
 279   int nlines = 0;
 280
 281   while ((getc (hmmfile)) != EOF) {
 282     int i;
 283
 284     if (nlines > MAXSEQLEN) {
 285       fprintf (stderr, "ERROR! the HMMER profile is too long. Can't cope... must quit.\n");
 286       exit (EXIT_FAILURE);
 287     }
 288
 289     for (i = 0; i < 24; i++) {
 290       fscanf (hmmfile, "%f", &data->hmmmat[nlines][i]);
 291     }
 292     nlines++;
 293   }
 294   data->lens = nlines - 1;      /* set the query sequence length */
 295 }
 296
 297 /*************************************************************************/
 298 /* Function to test whether the data array for input
 299    into the Neural Network function is of the expected size */
 300 void
 301 test_size (int size, int required_size, char *name)
 302 {
 303   if (size != required_size) {
 304     fprintf (stderr, "ERROR! input array not of correct size for %s (%i vs. %i)\n", name, size, required_size);
 305     exit (EXIT_FAILURE);
 306   }
 307 }
 308
 309
 310
 311 /*************************************************************************/
 312 /***********************************************************************/
 313 /* Functions below here haven't been checked */
 314 void
 315 dowinsspsi (float seq2str[MAXSEQLEN][3], int win, int curpos, float winarss[30][4], const alldata * data)
 316 {
 317   int i, j;
 318
 319   for (j = 0, i = (curpos - ((win - 1) / 2)); i <= (curpos + ((win - 1) / 2)); i++, j++) {
 320
 321     if (i >= 0 && i < data->lens) {
 322       /* CC - changesd to make it more similar to network training */
 323
 324       /* for (k = 0; k < 3; k++)
 325          winarss[j][k] = seq2str[i][k]; */
 326
 327       if ((seq2str[i][0] > seq2str[i][1]) && (seq2str[i][0] > seq2str[i][2])) {
 328         winarss[j][0] = 1.0;
 329         winarss[j][1] = 0.0;
 330         winarss[j][2] = 0.0;
 331       } else if ((seq2str[i][1] > seq2str[i][0]) && (seq2str[i][1] > seq2str[i][2])) {
 332         winarss[j][0] = 0.0;
 333         winarss[j][1] = 1.0;
 334         winarss[j][2] = 0.0;
 335       } else if ((seq2str[i][2] > seq2str[i][0]) && (seq2str[i][2] > seq2str[i][1])) {
 336         winarss[j][0] = 0.0;
 337         winarss[j][1] = 0.0;
 338         winarss[j][2] = 1.0;
 339       } else {
 340         winarss[j][0] = 0.0;
 341         winarss[j][1] = 0.0;
 342         winarss[j][2] = 1.0;
 343       }
 344     } else if (i < 0 || i >= data->lens) {
 345       /* Setting winarss[j][2] to 1 increases q3 accuracy */
 346       /* CC - does it really??? */
 347       winarss[j][0] = 0.0;
 348       winarss[j][1] = 0.0;
 349       winarss[j][2] = 0.0;
 350     }
 351   }
 352 }
 353
 354 /*************************************************************************/
 355 /* Structure to hold data from networks */
 356 typedef struct netdata
 357 {
 358   int ps1;
 359   float *netin;
 360   float *netprofin;
 361 } netdata;
 362
 363 /*************************************************************************/
 364 /* main prediction routine */
 365 void
 366 dopred (const alldata * data)
 367 {
 368   extern int NOHMM, NOPSI;
 369
 370   float psi2net[MAXSEQLEN][3];
 371   float hmmnet[MAXSEQLEN][3];
 372
 373   float finalout[MAXSEQLEN][3];
 374   int psi2fin[MAXSEQLEN];
 375   int hmmfin[MAXSEQLEN];
 376   int consfin[MAXSEQLEN];
 377
 378   char psi2let[MAXSEQLEN];
 379   char hmmlet[MAXSEQLEN];
 380   char finlet[MAXSEQLEN];
 381
 382   char sollet25[MAXSEQLEN];
 383   char sollet5[MAXSEQLEN];
 384   char sollet0[MAXSEQLEN];
 385
 386   float netprofin3[500];
 387   float confidence[MAXSEQLEN];  /* Array for holding the confidence value calculated by doconf() */
 388   float psiar[WINAR][WINAR];
 389   int winar2[WINAR][WINAR];
 390
 391   float solacc25[MAXSEQLEN][2];
 392   float solacc5[MAXSEQLEN][2];
 393   float solacc0[MAXSEQLEN][2];
 394
 395   float winarss[WINAR][4];
 396
 397   char letfilt[MAXSEQLEN];
 398
 399   int i, t;
 400   char jury[MAXSEQLEN];
 401
 402   float netprofin[MAXSEQLEN];
 403   memset (netprofin, 0x00, sizeof (float) * MAXSEQLEN);
 404   memset (winar2, 0x00, sizeof (winar2[0][0]) * WINAR * WINAR);
 405
 406   letfilt[0] = '\0';
 407
 408 #ifdef HMM
 409   {
 410     float seq2str[MAXSEQLEN][3];
 411     int i, windows = 17;
 412     for (i = 0; i < data->lens; i++) {
 413       float nn_output[3];
 414
 415       int y, j = 0;
 416       doprofhmm (data, windows, i, psiar);
 417
 418       for (y = 0; y < windows; y++) {
 419         int z;
 420         for (z = 0; z < PROLEN; z++)
 421           netprofin3[j++] = psiar[y][z];
 422       }
 423
 424       test_size (j, hmm1REC.NoOfInput, "hmm1");
 425
 426       hmm1 (netprofin3, nn_output, 0);
 427       seq2str[i][0] = nn_output[0];
 428       seq2str[i][1] = nn_output[1];
 429       seq2str[i][2] = nn_output[2];
 430     }
 431
 432     windows = 19;
 433     for (i = 0; i < data->lens; i++) {
 434       int y, j = 0;
 435       float nn_output[3];
 436       float nn_input[19 * 3];
 437
 438       dowinsspsi (seq2str, windows, i, winarss, data);
 439
 440       for (y = 0; y < windows; y++) {
 441         int z;
 442         for (z = 0; z < 3; z++)
 443           nn_input[j++] = winarss[y][z];
 444       }
 445
 446       test_size (j, hmm2REC.NoOfInput, "hmm2");
 447
 448       hmm2 (nn_input, nn_output, 0);
 449       hmmnet[i][0] = nn_output[0];
 450       hmmnet[i][1] = nn_output[1];
 451       hmmnet[i][2] = nn_output[2];
 452     }
 453   }
 454 #endif /* end ifdef HMM */
 455
 456   if (NOHMM == 0) {
 457     int i, windows = 17;
 458     for (i = 0; i < data->lens; i++) {
 459       float nn_output[3];
 460       int y, j = 0;
 461       doprofhmm (data, windows, i, psiar);
 462
 463       for (y = 0; y < windows; y++) {
 464         int z;
 465         for (z = 0; z < 24; z++)
 466           netprofin3[j++] = psiar[y][z];
 467       }
 468
 469       test_size (j, hmmsol25REC.NoOfInput, "hmmsol25");
 470
 471       hmmsol25 (netprofin3, nn_output, 0);
 472       solacc25[i][0] = nn_output[0];
 473       solacc25[i][1] = nn_output[1];
 474
 475       test_size (j, hmmsol5REC.NoOfInput, "hmmsol5");
 476
 477       hmmsol5 (netprofin3, nn_output, 0);
 478       solacc5[i][0] = nn_output[0];
 479       solacc5[i][1] = nn_output[1];
 480
 481       test_size (j, hmmsol0REC.NoOfInput, "hmmsol0");
 482
 483       hmmsol0 (netprofin3, nn_output, 0);
 484       solacc0[i][0] = nn_output[0];
 485       solacc0[i][1] = nn_output[1];
 486     }
 487   }
 488
 489   if (NOPSI == 0) {
 490     int i, windows = 17;
 491     for (i = 0; i < data->lens; i++) {
 492       float nn_output[3];
 493       int y, j = 0;
 494       doprofpsi (data, windows, i, psiar);
 495
 496       for (y = 0; y < windows; y++) {
 497         int z;
 498         for (z = 0; z < 20; z++)
 499           netprofin3[j++] = psiar[y][z];
 500       }
 501
 502       test_size (j, psisol25REC.NoOfInput, "psisol25");
 503
 504       psisol25 (netprofin3, nn_output, 0);
 505       solacc25[i][0] += nn_output[0];
 506       solacc25[i][1] += nn_output[1];
 507
 508       test_size (j, psisol5REC.NoOfInput, "psisol5");
 509
 510       psisol5 (netprofin3, nn_output, 0);
 511       solacc5[i][0] += nn_output[0];
 512       solacc5[i][1] += nn_output[1];
 513
 514       test_size (j, psisol0REC.NoOfInput, "psisol0");
 515
 516       psisol0 (netprofin3, nn_output, 0);
 517       solacc0[i][0] += nn_output[0];
 518       solacc0[i][1] += nn_output[1];
 519     }
 520   }
 521
 522   for (i = 0; i < data->lens; i++) {
 523     if (solacc25[i][0] > solacc25[i][1]) {
 524       sollet25[i] = '-';
 525     }
 526     if (solacc25[i][1] > solacc25[i][0]) {
 527       sollet25[i] = 'B';
 528     }
 529     if (solacc5[i][0] > solacc5[i][1]) {
 530       sollet5[i] = '-';
 531     }
 532     if (solacc5[i][1] > solacc5[i][0]) {
 533       sollet5[i] = 'B';
 534     }
 535     if (solacc0[i][0] > solacc0[i][1]) {
 536       sollet0[i] = '-';
 537     }
 538     if (solacc0[i][1] > solacc0[i][0]) {
 539       sollet0[i] = 'B';
 540     }
 541   }
 542
 543   if (NOPSI == 0) {
 544     float seq2str[MAXSEQLEN][3];
 545     int i, windows = 17;
 546
 547     /* This is for the PSIBLAST PSSM, for which there are two networks for each layer */
 548     memset (psi2net, 0x00, MAXSEQLEN * 3);
 549 #ifdef PSSM
 550     for (t = 0; t < 2; t++) {
 551
 552       /* First network */
 553       windows = 17;
 554       for (i = 0; i < data->lens; i++) {
 555         int y, j = 0;
 556         float nn_output[3];
 557
 558         doprofpsi (data, windows, i, psiar);
 559
 560         for (y = 0; y < windows; y++) {
 561           int z;
 562           for (z = 0; z < 20; z++)
 563             netprofin[j++] = psiar[y][z];
 564         }
 565
 566
 567         if (t == 0) {
 568           test_size (j, pssma1REC.NoOfInput, "pssma1");
 569           pssma1 (netprofin, nn_output, 0);
 570         } else if (t == 1) {
 571           test_size (j, pssmb1REC.NoOfInput, "pssmb1");
 572           pssmb1 (netprofin, nn_output, 0);
 573         }
 574
 575         seq2str[i][0] = nn_output[0];
 576         seq2str[i][1] = nn_output[1];
 577         seq2str[i][2] = nn_output[2];
 578       }
 579
 580       /* Second network */
 581       windows = 19;
 582       for (i = 0; i < data->lens; i++) {
 583         int y, j = 0;
 584         float nn_output[3];
 585         float nn_input[19 * 3];
 586
 587         dowinsspsi (seq2str, windows, i, winarss, data);
 588
 589         for (y = 0; y < windows; y++) {
 590           int z;
 591           for (z = 0; z < 3; z++)
 592             nn_input[j++] = winarss[y][z];
 593         }
 594
 595
 596         if (t == 0) {
 597           test_size (j, pssma2REC.NoOfInput, "pssma2");
 598           pssma2 (nn_input, nn_output, 0);
 599         } else if (t == 1) {
 600           test_size (j, pssmb2REC.NoOfInput, "pssmb2");
 601           pssmb2 (nn_input, nn_output, 0);
 602         }
 603
 604         psi2net[i][0] += nn_output[0];
 605         psi2net[i][1] += nn_output[1];
 606         psi2net[i][2] += nn_output[2];
 607       }
 608     }
 609     for (i = 0; i < data->lens; i++) {
 610       psi2net[i][0] /= 2;
 611       psi2net[i][1] /= 2;
 612       psi2net[i][2] /= 2;
 613     }
 614 #else
 615     for (i = 0; i < data->lens; i++) {
 616       psi2net[i][0] = 0;
 617       psi2net[i][1] = 0;
 618       psi2net[i][2] = 0;
 619     }
 620 #endif /* end of ifdef PSSM */
 621   }
 622
 623   /* Work out the final output values for the various outputs
 624    * There is some smoothing going on at the ends of the sequences
 625    * Also, some of these are being done even if there is no data in them.
 626    * */
 627
 628   if (NOPSI == 0) {
 629 #ifdef PSSM
 630     for (i = 0; i < data->lens; i++) {
 631       finalout[i][0] = psi2net[i][0];
 632       finalout[i][1] = psi2net[i][1];
 633       finalout[i][2] = psi2net[i][2];
 634
 635 #ifdef POST
 636       if (i <= 4)
 637         finalout[i][2] += (5 - i) * 0.2;
 638       else if (i >= (data->lens - 5))
 639         finalout[i][2] += (5 - ((data->lens - 1) - i)) * 0.2;
 640 #endif
 641
 642       if (finalout[i][0] > finalout[i][1] && finalout[i][0] > finalout[i][2])
 643         psi2fin[i] = SHEET;
 644       if (finalout[i][1] > finalout[i][0] && finalout[i][1] > finalout[i][2])
 645         psi2fin[i] = HELIX;
 646       if (finalout[i][2] > finalout[i][0] && finalout[i][2] > finalout[i][1])
 647         psi2fin[i] = COIL;
 648     }
 649     int2pred (psi2fin, psi2let, data->lens);
 650 #endif /* end PSSM */
 651   }
 652 #ifdef HMM
 653   if (NOHMM == 0) {
 654     for (i = 0; i < data->lens; i++) {
 655       finalout[i][0] = hmmnet[i][0];
 656       finalout[i][1] = hmmnet[i][1];
 657       finalout[i][2] = hmmnet[i][2];
 658
 659 #ifdef POST
 660       if (i <= 4)
 661         finalout[i][2] += (5 - i) * 0.2;
 662       else if (i >= (data->lens - 5))
 663         finalout[i][2] += (5 - ((data->lens - 1) - i)) * 0.2;
 664 #endif
 665
 666       if (finalout[i][0] > finalout[i][1] && finalout[i][0] > finalout[i][2])
 667         hmmfin[i] = SHEET;
 668       if (finalout[i][1] > finalout[i][0] && finalout[i][1] > finalout[i][2])
 669         hmmfin[i] = HELIX;
 670       if (finalout[i][2] > finalout[i][0] && finalout[i][2] > finalout[i][1])
 671         hmmfin[i] = COIL;
 672     }
 673     int2pred (hmmfin, hmmlet, data->lens);
 674   }
 675 #endif /* end ifdef HMM */
 676
 677   /* When all data is available use the 'cons' network to decide all non-jury positions.
 678    * Also calculate the confidence values */
 679   if (NOPSI == 0) {
 680     int j;
 681     for (j = 0; j < data->lens; j++) {
 682
 683       /* Is there disagreement? */
 684       if (psi2let[j] == hmmlet[j]) {    /* No - use the PSI networks */
 685         jury[j] = ' ';
 686         consfin[j] = psi2fin[j];
 687
 688         confidence[j] = doconf (psi2net[j][0], psi2net[j][1], psi2net[j][2]);
 689
 690       } else {                  /* yes - take the mean propensities of the HMM and PSI networks */
 691         jury[j] = '*';
 692
 693 #ifdef HMMONLY
 694         finalout[j][0] = hmmnet[j][0];
 695         finalout[j][1] = hmmnet[j][1];
 696         finalout[j][2] = hmmnet[j][2];
 697 #else
 698         finalout[j][0] = (hmmnet[j][0] + psi2net[j][0]) / 2;
 699         finalout[j][1] = (hmmnet[j][1] + psi2net[j][1]) / 2;
 700         finalout[j][2] = (hmmnet[j][2] + psi2net[j][2]) / 2;
 701 #endif /* end HMMONLY */
 702
 703         confidence[j] = doconf (finalout[j][0], finalout[j][1], finalout[j][2]);
 704
 705         consfin[j] = COIL;
 706
 707         if (finalout[j][0] > finalout[j][1] && finalout[j][0] > finalout[j][2])
 708           consfin[j] = SHEET;
 709         if (finalout[j][1] > finalout[j][0] && finalout[j][1] > finalout[j][2])
 710           consfin[j] = HELIX;
 711         if (finalout[j][2] > finalout[j][0] && finalout[j][2] > finalout[j][1])
 712           consfin[j] = COIL;
 713
 714       }                         /* agreement if/else block */
 715
 716     }                           /* end of sequence */
 717     consfin[data->lens] = 25;
 718   }
 719
 720   /* end of if (NOHMM==0 && NOPSI==0) */
 721   /* Calculate the confidence values */
 722   if (NOPSI == 1) {
 723     fprintf (stdout, "\n\nWARNING!: Only using the HMM profile\nAccuracy will average 79.6%%\n\n");
 724     for (i = 0; i < data->lens; i++) {
 725       consfin[i] = hmmfin[i];
 726
 727       confidence[i] = doconf ((hmmnet[i][0]), (hmmnet[i][1]), (hmmnet[i][2]));
 728     }
 729     consfin[data->lens] = 25;
 730   }
 731
 732   if (NOPSI == 0) {
 733     if (1 == TESTOUTPUT) {
 734       fprintf (stderr, "\n\nBoth PSIBLAST and HMM profiles were found\nAccuracy will average 81.6%%\n\n");
 735     } else {
 736       fprintf (stdout, "\n\nBoth PSIBLAST and HMM profiles were found\nAccuracy will average 81.6%%\n\n");
 737     }
 738   }
 739   int2pred (consfin, finlet, data->lens);
 740   finlet[data->lens] = '\0';
 741
 742 #ifdef FILTER
 743   /* Remove unlikely secondary structure */
 744   filter (finlet, letfilt);
 745   filter (finlet, letfilt);
 746   filter (finlet, letfilt);
 747 #else
 748   memcpy (letfilt, finlet, sizeof (char) * MAXSEQLEN);
 749 #endif
 750
 751   /* output predictions */
 752   if (HUMAN)                    /* print human readable and quit */
 753     print_human (data->lens, letfilt, confidence);
 754
 755   printf ("\njnetpred:");
 756   for (i = 0; i < data->lens; i++)
 757     printf ("%c,", letfilt[i]);
 758   printf ("\nJNETCONF:");
 759   for (i = 0; i < data->lens; i++)
 760     printf ("%1.0f,", confidence[i]);
 761   printf ("\nJNETSOL25:");
 762   for (i = 0; i < data->lens; i++)
 763     printf ("%c,", sollet25[i]);
 764   printf ("\nJNETSOL5:");
 765   for (i = 0; i < data->lens; i++)
 766     printf ("%c,", sollet5[i]);
 767   printf ("\nJNETSOL0:");
 768   for (i = 0; i < data->lens; i++)
 769     printf ("%c,", sollet0[i]);
 770   printf ("\nJNETHMM:");
 771   for (i = 0; i < data->lens; i++)
 772     printf ("%c,", hmmlet[i]);
 773
 774   if (NOPSI == 0) {
 775     printf ("\nJNETPSSM:");
 776     for (i = 0; i < data->lens; i++)
 777       printf ("%c,", psi2let[i]);
 778     printf ("\nJNETJURY:");
 779     for (i = 0; i < data->lens; i++)
 780       printf ("%c,", jury[i]);
 781   }
 782
 783   printf ("\nJNETPROPE:");
 784   for (i = 0; i < data->lens; i++)
 785     printf ("%.4f,", hmmnet[i][0]);
 786   printf ("\nJNETPROPH:");
 787   for (i = 0; i < data->lens; i++)
 788     printf ("%.4f,", hmmnet[i][1]);
 789   printf ("\nJNETPROPC:");
 790   for (i = 0; i < data->lens; i++)
 791     printf ("%.4f,", hmmnet[i][2]);
 792   printf ("\n");
 793 }
 794
 795
 796 /*************************************************************************/
 797 /* Calculate confidence as for PHD method (Rost and Sander, JMB, 1993, v232, p584) */
 798 float
 799 doconf (const float confa, const float confb, const float confc)
 800 {
 801   float whichout, maxout, maxnext, outconf;
 802   whichout = outconf = maxnext = 0;
 803   maxout = confc;
 804
 805   if (confa > confb && confa > confc) {
 806     whichout = 0;
 807     maxout = confa;
 808   }
 809   if (confb > confa && confb > confc) {
 810     whichout = 1;
 811     maxout = confb;
 812   }
 813   if (confc > confa && confc > confb) {
 814     whichout = 2;
 815     maxout = confc;
 816   }
 817
 818   if (whichout == 0) {
 819     if (confb > confc)
 820       maxnext = confb;
 821     if (confc > confb)
 822       maxnext = confc;
 823   }
 824   if (whichout == 1) {
 825     if (confc > confa)
 826       maxnext = confc;
 827     if (confa > confc)
 828       maxnext = confa;
 829   }
 830   if (whichout == 2) {
 831     if (confb > confa)
 832       maxnext = confb;
 833     if (confa > confb)
 834       maxnext = confa;
 835   }
 836   outconf = 10 * (maxout - maxnext);
 837   if (outconf > 9)
 838     outconf = 9;
 839
 840   return outconf;
 841 }
 842
 843 /*************************************************************************/
 844 void
 845 StrReplStr (char *s1, char *s2, char *FromSStr, char *ToSStr)
 846 {
 847   char *ChP1, *ChP2;
 848
 849   s1[0] = '\0';
 850   ChP1 = s2;
 851
 852   while ((ChP2 = strstr (ChP1, FromSStr)) != NULL) {
 853     if (ChP1 != ChP2)
 854       Strncat (s1, ChP1, strlen (ChP1) - strlen (ChP2));
 855     strcat (s1, ToSStr);
 856     ChP1 = ChP2 + strlen (FromSStr);
 857   }
 858   strcat (s1, ChP1);
 859   return;
 860 }
 861
 862 /*************************************************************************/
 863 char *
 864 Strncat (char *s1, char *s2, int len)
 865 {
 866   int OrigLen = 0;
 867   if (len == 0) {
 868     fprintf (stderr, "ERROR! Strncat error!");
 869     return s1;
 870   }
 871
 872   if (s1 == NULL || s2 == NULL) {
 873     fprintf (stderr, "ERROR! Strncat error!");
 874     return NULL;
 875   }
 876   OrigLen = strlen (s1);
 877   if (strncat (s1, s2, len) == NULL) {
 878     fprintf (stderr, "ERROR! Strncat error!");
 879     return NULL;
 880   }
 881
 882   s1[OrigLen + len] = '\0';
 883
 884   return s1;
 885 }
 886
 887 /*************************************************************************/
 888 /* Finds the data around the current position in a window */
 889 void
 890 doprofpsi (const alldata * data, const int win, const int curpos, float psiar[30][30])
 891 {
 892   int i, j = 0;
 893
 894   for (i = (curpos - ((win - 1) / 2)); i <= (curpos + ((win - 1) / 2)); i++) {
 895     int k;
 896     for (k = 0; k < 20; k++) {
 897       psiar[j][k] = data->psimat[i][k];
 898     }
 899
 900     if (i < 0 || i >= data->lens)
 901       memset (psiar[j], 0x00, sizeof (psiar[0][0]) * 30);
 902
 903     j++;
 904   }
 905 }
 906
 907 /*************************************************************************/
 908 /* Converts an array of predicted structure to a character array, not a null terminated string */
 909 void
 910 int2pred (const int *pred, char *letpred, const int length)
 911 {
 912   int i;
 913
 914   for (i = 0; i < length; i++) {
 915     switch (pred[i]) {
 916     case HELIX:
 917       letpred[i] = 'H';
 918       break;
 919     case SHEET:
 920       letpred[i] = 'E';
 921       break;
 922     case COIL:
 923       letpred[i] = '-';
 924       break;
 925       /*case 25: letpred[i] = '\0'; break; */
 926     default:
 927       fprintf (stderr, "ERROR! unknown secondary structure type '%i' at position %i in int2pred \n", pred[i], i);
 928       exit (EXIT_FAILURE);
 929       break;
 930     }
 931   }
 932 }
 933
 934 /*************************************************************************/
 935 void
 936 doprofhmm (const alldata * data, const int win, const int curpos, float psiar[30][30])
 937 {
 938   int i, j, k;
 939   j = k = i = 0;
 940
 941   for (i = (curpos - ((win - 1) / 2)); i <= (curpos + ((win - 1) / 2)); i++) {
 942     for (k = 0; k < 24; k++)
 943       psiar[j][k] = data->hmmmat[i][k];
 944
 945     if (i < 0 || i >= data->lens)
 946       for (k = 0; k < 24; k++)
 947         psiar[j][k] = 0;
 948     j++;
 949   }
 950 }
 951
 952 /*************************************************************************/
 953 /** print out a human readable version of the prediction and confidence **/
 954 void
 955 print_human (int length, char pred[MAXSEQLEN], float conf[MAXSEQLEN])
 956 {
 957   int i, j = 0, k;              /* counters */
 958
 959   float *confout;               /* array pointers to temp storage of output lines */
 960   char *predout;
 961
 962   /* allocate memory to array pointers */
 963   confout = malloc (BLOCK * sizeof (float));
 964   if (confout == NULL) {
 965     fprintf (stderr, "ERROR! Out of memory\n");
 966     exit (EXIT_FAILURE);
 967   }
 968
 969   predout = malloc (BLOCK * sizeof (char));
 970   if (predout == NULL) {
 971     fprintf (stderr, "ERROR! Out of memory\n");
 972     exit (EXIT_FAILURE);
 973   }
 974
 975   for (i = 0; i < length; i++) {
 976     confout[j] = conf[i];       /* store current conf value */
 977     predout[j] = pred[i];       /* store current prediction value */
 978
 979     /* when reached end of block, print stuff out */
 980     if (i > 1 && i % BLOCK == 0) {
 981       printf ("\n\npred:  ");
 982       for (k = 0; k < BLOCK; ++k)
 983         printf ("%c", predout[k]);
 984
 985       printf ("\nconf:  ");
 986       for (k = 0; k < BLOCK; ++k)
 987         printf ("%1.0f", confout[k]);
 988
 989       j = 0;                    /* zero the block counter */
 990     }
 991
 992     ++j;
 993   }
 994
 995   /* catch any trailing blocks that are smaller than required */
 996   if (j)
 997     printf ("\n\npred:  ");
 998   for (k = 0; k < j; ++k)
 999     printf ("%c", predout[k]);
1000
1001   printf ("\nconf:  ");
1002   for (k = 0; k < j; ++k)
1003     printf ("%1.0f", confout[k]);
1004   printf ("\n");
1005
1006   /* free memory and exit */
1007   free (confout);
1008   free (predout);
1009   exit (EXIT_SUCCESS);
1010 }