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     float tempval;
 260     int i;
 261     for (i = 0; i < 20; i++) {
 262       int npar = fscanf (psifile, "%f", &tempval);
 263       if (1 == npar) {
 264         data->psimat[nlines][i] = tempval;
 265       } else {
 266         fprintf (stderr, "ERROR: Can't read a value from PSSM file");
 267         data->psimat[nlines][i] = 0.;
 268       }
 269     }
 270     nlines++;
 271   }
 272
 273   /* check the data length matches that of the HMMer data */
 274   if (data->lens != nlines - 1) {
 275     fprintf (stderr, "ERROR! PSSM data length doesn't match that for HMMer. Are both input files from the same alignment/sequence?\n");
 276     exit (EXIT_FAILURE);
 277   }
 278 }
 279
 280
 281 /*************************************************************************/
 282 /* Reads in HMMER profile that was originally based upon ClustalW alignments */
 283 /* Additionally, use this data to define the length of the sequence in the absence of an alignment */
 284 void
 285 readhmm (FILE * hmmfile, alldata * data)
 286 {
 287   int nlines = 0;
 288
 289   while ((getc (hmmfile)) != EOF) {
 290     int i;
 291
 292     if (nlines > MAXSEQLEN) {
 293       fprintf (stderr, "ERROR! the HMMER profile is too long. Can't cope... must quit.\n");
 294       exit (EXIT_FAILURE);
 295     }
 296
 297     for (i = 0; i < 24; i++) {
 298       float tempval;
 299       int npar = fscanf (hmmfile, "%f", &tempval);
 300       if (1 == npar) {
 301         data->hmmmat[nlines][i] = tempval;
 302       } else {
 303         fprintf (stderr, "ERROR: Can't read a value from HMMER profile file");
 304         data->hmmmat[nlines][i] = 0.;
 305       }
 306     }
 307     nlines++;
 308   }
 309   data->lens = nlines - 1;      /* set the query sequence length */
 310 }
 311
 312 /*************************************************************************/
 313 /* Function to test whether the data array for input
 314    into the Neural Network function is of the expected size */
 315 void
 316 test_size (int size, int required_size, char *name)
 317 {
 318   if (size != required_size) {
 319     fprintf (stderr, "ERROR! input array not of correct size for %s (%i vs. %i)\n", name, size, required_size);
 320     exit (EXIT_FAILURE);
 321   }
 322 }
 323
 324
 325
 326 /*************************************************************************/
 327 /***********************************************************************/
 328 /* Functions below here haven't been checked */
 329 void
 330 dowinsspsi (float seq2str[MAXSEQLEN][3], int win, int curpos, float winarss[30][4], const alldata * data)
 331 {
 332   int i, j;
 333
 334   for (j = 0, i = (curpos - ((win - 1) / 2)); i <= (curpos + ((win - 1) / 2)); i++, j++) {
 335
 336     if (i >= 0 && i < data->lens) {
 337       /* CC - changesd to make it more similar to network training */
 338
 339       /* for (k = 0; k < 3; k++)
 340          winarss[j][k] = seq2str[i][k]; */
 341
 342       if ((seq2str[i][0] > seq2str[i][1]) && (seq2str[i][0] > seq2str[i][2])) {
 343         winarss[j][0] = 1.0;
 344         winarss[j][1] = 0.0;
 345         winarss[j][2] = 0.0;
 346       } else if ((seq2str[i][1] > seq2str[i][0]) && (seq2str[i][1] > seq2str[i][2])) {
 347         winarss[j][0] = 0.0;
 348         winarss[j][1] = 1.0;
 349         winarss[j][2] = 0.0;
 350       } else if ((seq2str[i][2] > seq2str[i][0]) && (seq2str[i][2] > seq2str[i][1])) {
 351         winarss[j][0] = 0.0;
 352         winarss[j][1] = 0.0;
 353         winarss[j][2] = 1.0;
 354       } else {
 355         winarss[j][0] = 0.0;
 356         winarss[j][1] = 0.0;
 357         winarss[j][2] = 1.0;
 358       }
 359     } else if (i < 0 || i >= data->lens) {
 360       /* Setting winarss[j][2] to 1 increases q3 accuracy */
 361       /* CC - does it really??? */
 362       winarss[j][0] = 0.0;
 363       winarss[j][1] = 0.0;
 364       winarss[j][2] = 0.0;
 365     }
 366   }
 367 }
 368
 369 /*************************************************************************/
 370 /* Structure to hold data from networks */
 371 typedef struct netdata
 372 {
 373   int ps1;
 374   float *netin;
 375   float *netprofin;
 376 } netdata;
 377
 378 /*************************************************************************/
 379 /* main prediction routine */
 380 void
 381 dopred (const alldata * data)
 382 {
 383   extern int NOHMM, NOPSI;
 384
 385   float psi2net[MAXSEQLEN][3];
 386   float hmmnet[MAXSEQLEN][3];
 387
 388   float finalout[MAXSEQLEN][3];
 389   int psi2fin[MAXSEQLEN];
 390   int hmmfin[MAXSEQLEN];
 391   int consfin[MAXSEQLEN];
 392
 393   char psi2let[MAXSEQLEN];
 394   char hmmlet[MAXSEQLEN];
 395   char finlet[MAXSEQLEN];
 396
 397   char sollet25[MAXSEQLEN];
 398   char sollet5[MAXSEQLEN];
 399   char sollet0[MAXSEQLEN];
 400
 401   float netprofin3[500];
 402   float confidence[MAXSEQLEN];  /* Array for holding the confidence value calculated by doconf() */
 403   float psiar[WINAR][WINAR];
 404   int winar2[WINAR][WINAR];
 405
 406   float solacc25[MAXSEQLEN][2];
 407   float solacc5[MAXSEQLEN][2];
 408   float solacc0[MAXSEQLEN][2];
 409
 410   float winarss[WINAR][4];
 411
 412   char letfilt[MAXSEQLEN];
 413
 414   int i, t;
 415   char jury[MAXSEQLEN];
 416
 417   float netprofin[MAXSEQLEN];
 418   memset (netprofin, 0x00, sizeof (float) * MAXSEQLEN);
 419   memset (winar2, 0x00, sizeof (winar2[0][0]) * WINAR * WINAR);
 420
 421   letfilt[0] = '\0';
 422
 423 #ifdef HMM
 424   {
 425     float seq2str[MAXSEQLEN][3];
 426     int i, windows = 17;
 427     for (i = 0; i < data->lens; i++) {
 428       float nn_output[3];
 429
 430       int y, j = 0;
 431       doprofhmm (data, windows, i, psiar);
 432
 433       for (y = 0; y < windows; y++) {
 434         int z;
 435         for (z = 0; z < PROLEN; z++)
 436           netprofin3[j++] = psiar[y][z];
 437       }
 438
 439       test_size (j, hmm1REC.NoOfInput, "hmm1");
 440
 441       hmm1 (netprofin3, nn_output, 0);
 442       seq2str[i][0] = nn_output[0];
 443       seq2str[i][1] = nn_output[1];
 444       seq2str[i][2] = nn_output[2];
 445     }
 446
 447     windows = 19;
 448     for (i = 0; i < data->lens; i++) {
 449       int y, j = 0;
 450       float nn_output[3];
 451       float nn_input[19 * 3];
 452
 453       dowinsspsi (seq2str, windows, i, winarss, data);
 454
 455       for (y = 0; y < windows; y++) {
 456         int z;
 457         for (z = 0; z < 3; z++)
 458           nn_input[j++] = winarss[y][z];
 459       }
 460
 461       test_size (j, hmm2REC.NoOfInput, "hmm2");
 462
 463       hmm2 (nn_input, nn_output, 0);
 464       hmmnet[i][0] = nn_output[0];
 465       hmmnet[i][1] = nn_output[1];
 466       hmmnet[i][2] = nn_output[2];
 467     }
 468   }
 469 #endif /* end ifdef HMM */
 470
 471   if (NOHMM == 0) {
 472     int i, windows = 17;
 473     for (i = 0; i < data->lens; i++) {
 474       float nn_output[3];
 475       int y, j = 0;
 476       doprofhmm (data, windows, i, psiar);
 477
 478       for (y = 0; y < windows; y++) {
 479         int z;
 480         for (z = 0; z < 24; z++)
 481           netprofin3[j++] = psiar[y][z];
 482       }
 483
 484       test_size (j, hmmsol25REC.NoOfInput, "hmmsol25");
 485
 486       hmmsol25 (netprofin3, nn_output, 0);
 487       solacc25[i][0] = nn_output[0];
 488       solacc25[i][1] = nn_output[1];
 489
 490       test_size (j, hmmsol5REC.NoOfInput, "hmmsol5");
 491
 492       hmmsol5 (netprofin3, nn_output, 0);
 493       solacc5[i][0] = nn_output[0];
 494       solacc5[i][1] = nn_output[1];
 495
 496       test_size (j, hmmsol0REC.NoOfInput, "hmmsol0");
 497
 498       hmmsol0 (netprofin3, nn_output, 0);
 499       solacc0[i][0] = nn_output[0];
 500       solacc0[i][1] = nn_output[1];
 501     }
 502   }
 503
 504   if (NOPSI == 0) {
 505     int i, windows = 17;
 506     for (i = 0; i < data->lens; i++) {
 507       float nn_output[3];
 508       int y, j = 0;
 509       doprofpsi (data, windows, i, psiar);
 510
 511       for (y = 0; y < windows; y++) {
 512         int z;
 513         for (z = 0; z < 20; z++)
 514           netprofin3[j++] = psiar[y][z];
 515       }
 516
 517       test_size (j, psisol25REC.NoOfInput, "psisol25");
 518
 519       psisol25 (netprofin3, nn_output, 0);
 520       solacc25[i][0] += nn_output[0];
 521       solacc25[i][1] += nn_output[1];
 522
 523       test_size (j, psisol5REC.NoOfInput, "psisol5");
 524
 525       psisol5 (netprofin3, nn_output, 0);
 526       solacc5[i][0] += nn_output[0];
 527       solacc5[i][1] += nn_output[1];
 528
 529       test_size (j, psisol0REC.NoOfInput, "psisol0");
 530
 531       psisol0 (netprofin3, nn_output, 0);
 532       solacc0[i][0] += nn_output[0];
 533       solacc0[i][1] += nn_output[1];
 534     }
 535   }
 536
 537   for (i = 0; i < data->lens; i++) {
 538     if (solacc25[i][0] > solacc25[i][1]) {
 539       sollet25[i] = '-';
 540     }
 541     if (solacc25[i][1] > solacc25[i][0]) {
 542       sollet25[i] = 'B';
 543     }
 544     if (solacc5[i][0] > solacc5[i][1]) {
 545       sollet5[i] = '-';
 546     }
 547     if (solacc5[i][1] > solacc5[i][0]) {
 548       sollet5[i] = 'B';
 549     }
 550     if (solacc0[i][0] > solacc0[i][1]) {
 551       sollet0[i] = '-';
 552     }
 553     if (solacc0[i][1] > solacc0[i][0]) {
 554       sollet0[i] = 'B';
 555     }
 556   }
 557
 558   if (NOPSI == 0) {
 559     float seq2str[MAXSEQLEN][3];
 560     int i, windows = 17;
 561
 562     /* This is for the PSIBLAST PSSM, for which there are two networks for each layer */
 563     memset (psi2net, 0x00, MAXSEQLEN * 3);
 564 #ifdef PSSM
 565     for (t = 0; t < 2; t++) {
 566
 567       /* First network */
 568       windows = 17;
 569       for (i = 0; i < data->lens; i++) {
 570         int y, j = 0;
 571         float nn_output[3];
 572
 573         doprofpsi (data, windows, i, psiar);
 574
 575         for (y = 0; y < windows; y++) {
 576           int z;
 577           for (z = 0; z < 20; z++)
 578             netprofin[j++] = psiar[y][z];
 579         }
 580
 581
 582         if (t == 0) {
 583           test_size (j, pssma1REC.NoOfInput, "pssma1");
 584           pssma1 (netprofin, nn_output, 0);
 585         } else if (t == 1) {
 586           test_size (j, pssmb1REC.NoOfInput, "pssmb1");
 587           pssmb1 (netprofin, nn_output, 0);
 588         }
 589
 590         seq2str[i][0] = nn_output[0];
 591         seq2str[i][1] = nn_output[1];
 592         seq2str[i][2] = nn_output[2];
 593       }
 594
 595       /* Second network */
 596       windows = 19;
 597       for (i = 0; i < data->lens; i++) {
 598         int y, j = 0;
 599         float nn_output[3];
 600         float nn_input[19 * 3];
 601
 602         dowinsspsi (seq2str, windows, i, winarss, data);
 603
 604         for (y = 0; y < windows; y++) {
 605           int z;
 606           for (z = 0; z < 3; z++)
 607             nn_input[j++] = winarss[y][z];
 608         }
 609
 610
 611         if (t == 0) {
 612           test_size (j, pssma2REC.NoOfInput, "pssma2");
 613           pssma2 (nn_input, nn_output, 0);
 614         } else if (t == 1) {
 615           test_size (j, pssmb2REC.NoOfInput, "pssmb2");
 616           pssmb2 (nn_input, nn_output, 0);
 617         }
 618
 619         psi2net[i][0] += nn_output[0];
 620         psi2net[i][1] += nn_output[1];
 621         psi2net[i][2] += nn_output[2];
 622       }
 623     }
 624     for (i = 0; i < data->lens; i++) {
 625       psi2net[i][0] /= 2;
 626       psi2net[i][1] /= 2;
 627       psi2net[i][2] /= 2;
 628     }
 629 #else
 630     for (i = 0; i < data->lens; i++) {
 631       psi2net[i][0] = 0;
 632       psi2net[i][1] = 0;
 633       psi2net[i][2] = 0;
 634     }
 635 #endif /* end of ifdef PSSM */
 636   }
 637
 638   /* Work out the final output values for the various outputs
 639    * There is some smoothing going on at the ends of the sequences
 640    * Also, some of these are being done even if there is no data in them.
 641    * */
 642
 643   if (NOPSI == 0) {
 644 #ifdef PSSM
 645     for (i = 0; i < data->lens; i++) {
 646       finalout[i][0] = psi2net[i][0];
 647       finalout[i][1] = psi2net[i][1];
 648       finalout[i][2] = psi2net[i][2];
 649
 650 #ifdef POST
 651       if (i <= 4)
 652         finalout[i][2] += (5 - i) * 0.2;
 653       else if (i >= (data->lens - 5))
 654         finalout[i][2] += (5 - ((data->lens - 1) - i)) * 0.2;
 655 #endif
 656
 657       if (finalout[i][0] > finalout[i][1] && finalout[i][0] > finalout[i][2])
 658         psi2fin[i] = SHEET;
 659       if (finalout[i][1] > finalout[i][0] && finalout[i][1] > finalout[i][2])
 660         psi2fin[i] = HELIX;
 661       if (finalout[i][2] > finalout[i][0] && finalout[i][2] > finalout[i][1])
 662         psi2fin[i] = COIL;
 663     }
 664     int2pred (psi2fin, psi2let, data->lens);
 665 #endif /* end PSSM */
 666   }
 667 #ifdef HMM
 668   if (NOHMM == 0) {
 669     for (i = 0; i < data->lens; i++) {
 670       finalout[i][0] = hmmnet[i][0];
 671       finalout[i][1] = hmmnet[i][1];
 672       finalout[i][2] = hmmnet[i][2];
 673
 674 #ifdef POST
 675       if (i <= 4)
 676         finalout[i][2] += (5 - i) * 0.2;
 677       else if (i >= (data->lens - 5))
 678         finalout[i][2] += (5 - ((data->lens - 1) - i)) * 0.2;
 679 #endif
 680
 681       if (finalout[i][0] > finalout[i][1] && finalout[i][0] > finalout[i][2])
 682         hmmfin[i] = SHEET;
 683       if (finalout[i][1] > finalout[i][0] && finalout[i][1] > finalout[i][2])
 684         hmmfin[i] = HELIX;
 685       if (finalout[i][2] > finalout[i][0] && finalout[i][2] > finalout[i][1])
 686         hmmfin[i] = COIL;
 687     }
 688     int2pred (hmmfin, hmmlet, data->lens);
 689   }
 690 #endif /* end ifdef HMM */
 691
 692   /* When all data is available use the 'cons' network to decide all non-jury positions.
 693    * Also calculate the confidence values */
 694   if (NOPSI == 0) {
 695     int j;
 696     for (j = 0; j < data->lens; j++) {
 697
 698       /* Is there disagreement? */
 699       if (psi2let[j] == hmmlet[j]) {    /* No - use the PSI networks */
 700         jury[j] = ' ';
 701         consfin[j] = psi2fin[j];
 702
 703         confidence[j] = doconf (psi2net[j][0], psi2net[j][1], psi2net[j][2]);
 704
 705       } else {                  /* yes - take the mean propensities of the HMM and PSI networks */
 706         jury[j] = '*';
 707
 708 #ifdef HMMONLY
 709         finalout[j][0] = hmmnet[j][0];
 710         finalout[j][1] = hmmnet[j][1];
 711         finalout[j][2] = hmmnet[j][2];
 712 #else
 713         finalout[j][0] = (hmmnet[j][0] + psi2net[j][0]) / 2;
 714         finalout[j][1] = (hmmnet[j][1] + psi2net[j][1]) / 2;
 715         finalout[j][2] = (hmmnet[j][2] + psi2net[j][2]) / 2;
 716 #endif /* end HMMONLY */
 717
 718         confidence[j] = doconf (finalout[j][0], finalout[j][1], finalout[j][2]);
 719
 720         consfin[j] = COIL;
 721
 722         if (finalout[j][0] > finalout[j][1] && finalout[j][0] > finalout[j][2])
 723           consfin[j] = SHEET;
 724         if (finalout[j][1] > finalout[j][0] && finalout[j][1] > finalout[j][2])
 725           consfin[j] = HELIX;
 726         if (finalout[j][2] > finalout[j][0] && finalout[j][2] > finalout[j][1])
 727           consfin[j] = COIL;
 728
 729       }                         /* agreement if/else block */
 730
 731     }                           /* end of sequence */
 732     consfin[data->lens] = 25;
 733   }
 734
 735   /* end of if (NOHMM==0 && NOPSI==0) */
 736   /* Calculate the confidence values */
 737   if (NOPSI == 1) {
 738     fprintf (stdout, "\n\nWARNING!: Only using the HMM profile\nAccuracy will average 79.6%%\n\n");
 739     for (i = 0; i < data->lens; i++) {
 740       consfin[i] = hmmfin[i];
 741
 742       confidence[i] = doconf ((hmmnet[i][0]), (hmmnet[i][1]), (hmmnet[i][2]));
 743     }
 744     consfin[data->lens] = 25;
 745   }
 746
 747   if (NOPSI == 0) {
 748     if (1 == TESTOUTPUT) {
 749       fprintf (stderr, "\n\nBoth PSIBLAST and HMM profiles were found\nAccuracy will average 81.6%%\n\n");
 750     } else {
 751       fprintf (stdout, "\n\nBoth PSIBLAST and HMM profiles were found\nAccuracy will average 81.6%%\n\n");
 752     }
 753   }
 754   int2pred (consfin, finlet, data->lens);
 755   finlet[data->lens] = '\0';
 756
 757 #ifdef FILTER
 758   /* Remove unlikely secondary structure */
 759   filter (finlet, letfilt);
 760   filter (finlet, letfilt);
 761   filter (finlet, letfilt);
 762 #else
 763   memcpy (letfilt, finlet, sizeof (char) * MAXSEQLEN);
 764 #endif
 765
 766   /* output predictions */
 767   if (HUMAN)                    /* print human readable and quit */
 768     print_human (data->lens, letfilt, confidence);
 769
 770   printf ("\njnetpred:");
 771   for (i = 0; i < data->lens; i++)
 772     printf ("%c,", letfilt[i]);
 773   printf ("\nJNETCONF:");
 774   for (i = 0; i < data->lens; i++)
 775     printf ("%1.0f,", confidence[i]);
 776   printf ("\nJNETSOL25:");
 777   for (i = 0; i < data->lens; i++)
 778     printf ("%c,", sollet25[i]);
 779   printf ("\nJNETSOL5:");
 780   for (i = 0; i < data->lens; i++)
 781     printf ("%c,", sollet5[i]);
 782   printf ("\nJNETSOL0:");
 783   for (i = 0; i < data->lens; i++)
 784     printf ("%c,", sollet0[i]);
 785   printf ("\nJNETHMM:");
 786   for (i = 0; i < data->lens; i++)
 787     printf ("%c,", hmmlet[i]);
 788
 789   if (NOPSI == 0) {
 790     printf ("\nJNETPSSM:");
 791     for (i = 0; i < data->lens; i++)
 792       printf ("%c,", psi2let[i]);
 793     printf ("\nJNETJURY:");
 794     for (i = 0; i < data->lens; i++)
 795       printf ("%c,", jury[i]);
 796   }
 797
 798   printf ("\nJNETPROPE:");
 799   for (i = 0; i < data->lens; i++)
 800     printf ("%.4f,", hmmnet[i][0]);
 801   printf ("\nJNETPROPH:");
 802   for (i = 0; i < data->lens; i++)
 803     printf ("%.4f,", hmmnet[i][1]);
 804   printf ("\nJNETPROPC:");
 805   for (i = 0; i < data->lens; i++)
 806     printf ("%.4f,", hmmnet[i][2]);
 807   printf ("\n");
 808 }
 809
 810
 811 /*************************************************************************/
 812 /* Calculate confidence as for PHD method (Rost and Sander, JMB, 1993, v232, p584) */
 813 float
 814 doconf (const float confa, const float confb, const float confc)
 815 {
 816   float whichout, maxout, maxnext, outconf;
 817   whichout = outconf = maxnext = 0;
 818   maxout = confc;
 819
 820   if (confa > confb && confa > confc) {
 821     whichout = 0;
 822     maxout = confa;
 823   }
 824   if (confb > confa && confb > confc) {
 825     whichout = 1;
 826     maxout = confb;
 827   }
 828   if (confc > confa && confc > confb) {
 829     whichout = 2;
 830     maxout = confc;
 831   }
 832
 833   if (whichout == 0) {
 834     if (confb > confc)
 835       maxnext = confb;
 836     if (confc > confb)
 837       maxnext = confc;
 838   }
 839   if (whichout == 1) {
 840     if (confc > confa)
 841       maxnext = confc;
 842     if (confa > confc)
 843       maxnext = confa;
 844   }
 845   if (whichout == 2) {
 846     if (confb > confa)
 847       maxnext = confb;
 848     if (confa > confb)
 849       maxnext = confa;
 850   }
 851   outconf = 10 * (maxout - maxnext);
 852   if (outconf > 9)
 853     outconf = 9;
 854
 855   return outconf;
 856 }
 857
 858 /*************************************************************************/
 859 void
 860 StrReplStr (char *s1, char *s2, char *FromSStr, char *ToSStr)
 861 {
 862   char *ChP1, *ChP2;
 863
 864   s1[0] = '\0';
 865   ChP1 = s2;
 866
 867   while ((ChP2 = strstr (ChP1, FromSStr)) != NULL) {
 868     if (ChP1 != ChP2)
 869       Strncat (s1, ChP1, strlen (ChP1) - strlen (ChP2));
 870     strcat (s1, ToSStr);
 871     ChP1 = ChP2 + strlen (FromSStr);
 872   }
 873   strcat (s1, ChP1);
 874   return;
 875 }
 876
 877 /*************************************************************************/
 878 char *
 879 Strncat (char *s1, char *s2, int len)
 880 {
 881   int OrigLen = 0;
 882   if (len == 0) {
 883     fprintf (stderr, "ERROR! Strncat error!");
 884     return s1;
 885   }
 886
 887   if (s1 == NULL || s2 == NULL) {
 888     fprintf (stderr, "ERROR! Strncat error!");
 889     return NULL;
 890   }
 891   OrigLen = strlen (s1);
 892   if (strncat (s1, s2, len) == NULL) {
 893     fprintf (stderr, "ERROR! Strncat error!");
 894     return NULL;
 895   }
 896
 897   s1[OrigLen + len] = '\0';
 898
 899   return s1;
 900 }
 901
 902 /*************************************************************************/
 903 /* Finds the data around the current position in a window */
 904 void
 905 doprofpsi (const alldata * data, const int win, const int curpos, float psiar[30][30])
 906 {
 907   int i, j = 0;
 908
 909   for (i = (curpos - ((win - 1) / 2)); i <= (curpos + ((win - 1) / 2)); i++) {
 910     int k;
 911     for (k = 0; k < 20; k++) {
 912       psiar[j][k] = data->psimat[i][k];
 913     }
 914
 915     if (i < 0 || i >= data->lens)
 916       memset (psiar[j], 0x00, sizeof (psiar[0][0]) * 30);
 917
 918     j++;
 919   }
 920 }
 921
 922 /*************************************************************************/
 923 /* Converts an array of predicted structure to a character array, not a null terminated string */
 924 void
 925 int2pred (const int *pred, char *letpred, const int length)
 926 {
 927   int i;
 928
 929   for (i = 0; i < length; i++) {
 930     switch (pred[i]) {
 931     case HELIX:
 932       letpred[i] = 'H';
 933       break;
 934     case SHEET:
 935       letpred[i] = 'E';
 936       break;
 937     case COIL:
 938       letpred[i] = '-';
 939       break;
 940       /*case 25: letpred[i] = '\0'; break; */
 941     default:
 942       fprintf (stderr, "ERROR! unknown secondary structure type '%i' at position %i in int2pred \n", pred[i], i);
 943       exit (EXIT_FAILURE);
 944       break;
 945     }
 946   }
 947 }
 948
 949 /*************************************************************************/
 950 void
 951 doprofhmm (const alldata * data, const int win, const int curpos, float psiar[30][30])
 952 {
 953   int i, j, k;
 954   j = k = i = 0;
 955
 956   for (i = (curpos - ((win - 1) / 2)); i <= (curpos + ((win - 1) / 2)); i++) {
 957     for (k = 0; k < 24; k++)
 958       psiar[j][k] = data->hmmmat[i][k];
 959
 960     if (i < 0 || i >= data->lens)
 961       for (k = 0; k < 24; k++)
 962         psiar[j][k] = 0;
 963     j++;
 964   }
 965 }
 966
 967 /*************************************************************************/
 968 /* print out a human readable version of the prediction and confidence */
 969 void
 970 print_human (int length, char pred[MAXSEQLEN], float conf[MAXSEQLEN])
 971 {
 972   int i, j = 0, k;
 973   float *confout;               /* array pointers to temp storage of output lines */
 974   char *predout;
 975
 976   /* allocate memory to array pointers */
 977   confout = malloc (BLOCK * sizeof (float));
 978   if (confout == NULL) {
 979     fprintf (stderr, "ERROR! Out of memory\n");
 980     exit (EXIT_FAILURE);
 981   }
 982
 983   predout = malloc (BLOCK * sizeof (char));
 984   if (predout == NULL) {
 985     fprintf (stderr, "ERROR! Out of memory\n");
 986     exit (EXIT_FAILURE);
 987   }
 988
 989   for (i = 0; i < length; i++) {
 990     confout[j] = conf[i];       /* store current conf value */
 991     predout[j] = pred[i];       /* store current prediction value */
 992
 993     /* when reached end of block, print stuff out */
 994     if (i > 1 && i % BLOCK == 0) {
 995       printf ("\n\npred:  ");
 996       for (k = 0; k < BLOCK; ++k)
 997         printf ("%c", predout[k]);
 998
 999       printf ("\nconf:  ");
1000       for (k = 0; k < BLOCK; ++k)
1001         printf ("%1.0f", confout[k]);
1002
1003       j = 0;                    /* zero the block counter */
1004     }
1005
1006     ++j;
1007   }
1008
1009   /* catch any trailing blocks that are smaller than required */
1010   if (j)
1011     printf ("\n\npred:  ");
1012   for (k = 0; k < j; ++k)
1013     printf ("%c", predout[k]);
1014
1015   printf ("\nconf:  ");
1016   for (k = 0; k < j; ++k)
1017     printf ("%1.0f", confout[k]);
1018   printf ("\n");
1019
1020   /* free memory and exit */
1021   free (confout);
1022   free (predout);
1023   exit (EXIT_SUCCESS);
1024 }