1 /************************************************************************
2 * Jnet - A consensus neural network protein *
3 * secondary structure prediction method *
5 * Copyright 1999,2009 James Cuff, Jonathan Barber and Christian Cole *
7 ************************************************************************
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.
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.
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 -------------------------------------------------------------------------*/
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
30 /*#define POST *//* Allow post processing of NN output */
31 /*#define FILTER *//* Allow string post processing of final output */
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 */
39 /* The allowed value of the integers representing residues */
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 */
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 */
60 #include "cmdline/cmdline.h"
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 */
67 #include "hmm1.h" /* First network for HMMER profile prediction */
68 #include "hmm2.h" /* Second network for hmm1 */
78 /* Main data structure for information about input data. Not all the fields are
80 typedef struct alldata
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 */
89 /* Global Vars - used to specify which input profiles have been provided */
93 int HUMAN = 0; /* toggle human readable output */
94 int BLOCK = 60; /* block size for human output */
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]);
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]);
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]);
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);
118 /*************************************************************/
120 main (int argc, char **argv)
122 FILE *psifile, *hmmfile;
125 /* Declared in cmdline/cmdline.h */
126 struct gengetopt_args_info args_info;
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__);
134 /* Get the command line arguments */
135 if (cmdline_parser (argc, argv, &args_info) != 0)
138 /* check whether TEST STYLE OUTPUT has been requested */
139 if (args_info.test_flag) {
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);
147 fprintf (stdout, "%s %s\n", CMDLINE_PARSER_PACKAGE, CMDLINE_PARSER_VERSION);
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);
155 if (1 == TESTOUTPUT) {
156 fprintf (stderr, "Found HMM profile data\n");
158 fprintf (stdout, "Found HMM profile data\n");
160 readhmm (hmmfile, data);
161 (void) fclose (hmmfile);
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);
170 if (1 == TESTOUTPUT) {
171 fprintf (stderr, "Found PSSM profile file\n");
173 fprintf (stdout, "Found PSSM profile file\n");
175 readpsi (psifile, data);
176 (void) fclose (psifile);
182 /* check whether human readable output is required (set by --human) */
183 if (args_info.human_given) {
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");
196 /* Do the prediction */
197 if (1 == TESTOUTPUT) {
198 fprintf (stderr, "Running final predictions!\n");
200 fprintf (stdout, "Running final predictions!\n");
210 /*************************************************************************/
211 /* Replaces particular runs of secondary structure with something else */
213 filter (char inp[MAXSEQLEN], char outp[MAXSEQLEN])
245 for (i = 0; pairs[i].orig != NULL; i++) {
246 StrReplStr (outp, inp, pairs[i].orig, pairs[i].final);
251 /*************************************************************************/
252 /* Reads in the PSIBLAST PSSM */
254 readpsi (FILE * psifile, alldata * data)
258 while ((getc (psifile)) != EOF) {
260 for (i = 0; i < 20; i++)
261 fscanf (psifile, "%f", &data->psimat[nlines][i]);
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");
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 */
277 readhmm (FILE * hmmfile, alldata * data)
281 while ((getc (hmmfile)) != EOF) {
284 if (nlines > MAXSEQLEN) {
285 fprintf (stderr, "ERROR! the HMMER profile is too long. Can't cope... must quit.\n");
289 for (i = 0; i < 24; i++) {
290 fscanf (hmmfile, "%f", &data->hmmmat[nlines][i]);
294 data->lens = nlines - 1; /* set the query sequence length */
297 /*************************************************************************/
298 /* Function to test whether the data array for input
299 into the Neural Network function is of the expected size */
301 test_size (int size, int required_size, char *name)
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);
311 /*************************************************************************/
312 /***********************************************************************/
313 /* Functions below here haven't been checked */
315 dowinsspsi (float seq2str[MAXSEQLEN][3], int win, int curpos, float winarss[30][4], const alldata * data)
319 for (j = 0, i = (curpos - ((win - 1) / 2)); i <= (curpos + ((win - 1) / 2)); i++, j++) {
321 if (i >= 0 && i < data->lens) {
322 /* CC - changesd to make it more similar to network training */
324 /* for (k = 0; k < 3; k++)
325 winarss[j][k] = seq2str[i][k]; */
327 if ((seq2str[i][0] > seq2str[i][1]) && (seq2str[i][0] > seq2str[i][2])) {
331 } else if ((seq2str[i][1] > seq2str[i][0]) && (seq2str[i][1] > seq2str[i][2])) {
335 } else if ((seq2str[i][2] > seq2str[i][0]) && (seq2str[i][2] > seq2str[i][1])) {
344 } else if (i < 0 || i >= data->lens) {
345 /* Setting winarss[j][2] to 1 increases q3 accuracy */
346 /* CC - does it really??? */
354 /*************************************************************************/
355 /* Structure to hold data from networks */
356 typedef struct netdata
363 /*************************************************************************/
364 /* main prediction routine */
366 dopred (const alldata * data)
368 extern int NOHMM, NOPSI;
370 float psi2net[MAXSEQLEN][3];
371 float hmmnet[MAXSEQLEN][3];
373 float finalout[MAXSEQLEN][3];
374 int psi2fin[MAXSEQLEN];
375 int hmmfin[MAXSEQLEN];
376 int consfin[MAXSEQLEN];
378 char psi2let[MAXSEQLEN];
379 char hmmlet[MAXSEQLEN];
380 char finlet[MAXSEQLEN];
382 char sollet25[MAXSEQLEN];
383 char sollet5[MAXSEQLEN];
384 char sollet0[MAXSEQLEN];
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];
391 float solacc25[MAXSEQLEN][2];
392 float solacc5[MAXSEQLEN][2];
393 float solacc0[MAXSEQLEN][2];
395 float winarss[WINAR][4];
397 char letfilt[MAXSEQLEN];
400 char jury[MAXSEQLEN];
402 float netprofin[MAXSEQLEN];
403 memset (netprofin, 0x00, sizeof (float) * MAXSEQLEN);
404 memset (winar2, 0x00, sizeof (winar2[0][0]) * WINAR * WINAR);
410 float seq2str[MAXSEQLEN][3];
412 for (i = 0; i < data->lens; i++) {
416 doprofhmm (data, windows, i, psiar);
418 for (y = 0; y < windows; y++) {
420 for (z = 0; z < PROLEN; z++)
421 netprofin3[j++] = psiar[y][z];
424 test_size (j, hmm1REC.NoOfInput, "hmm1");
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];
433 for (i = 0; i < data->lens; i++) {
436 float nn_input[19 * 3];
438 dowinsspsi (seq2str, windows, i, winarss, data);
440 for (y = 0; y < windows; y++) {
442 for (z = 0; z < 3; z++)
443 nn_input[j++] = winarss[y][z];
446 test_size (j, hmm2REC.NoOfInput, "hmm2");
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];
454 #endif /* end ifdef HMM */
458 for (i = 0; i < data->lens; i++) {
461 doprofhmm (data, windows, i, psiar);
463 for (y = 0; y < windows; y++) {
465 for (z = 0; z < 24; z++)
466 netprofin3[j++] = psiar[y][z];
469 test_size (j, hmmsol25REC.NoOfInput, "hmmsol25");
471 hmmsol25 (netprofin3, nn_output, 0);
472 solacc25[i][0] = nn_output[0];
473 solacc25[i][1] = nn_output[1];
475 test_size (j, hmmsol5REC.NoOfInput, "hmmsol5");
477 hmmsol5 (netprofin3, nn_output, 0);
478 solacc5[i][0] = nn_output[0];
479 solacc5[i][1] = nn_output[1];
481 test_size (j, hmmsol0REC.NoOfInput, "hmmsol0");
483 hmmsol0 (netprofin3, nn_output, 0);
484 solacc0[i][0] = nn_output[0];
485 solacc0[i][1] = nn_output[1];
491 for (i = 0; i < data->lens; i++) {
494 doprofpsi (data, windows, i, psiar);
496 for (y = 0; y < windows; y++) {
498 for (z = 0; z < 20; z++)
499 netprofin3[j++] = psiar[y][z];
502 test_size (j, psisol25REC.NoOfInput, "psisol25");
504 psisol25 (netprofin3, nn_output, 0);
505 solacc25[i][0] += nn_output[0];
506 solacc25[i][1] += nn_output[1];
508 test_size (j, psisol5REC.NoOfInput, "psisol5");
510 psisol5 (netprofin3, nn_output, 0);
511 solacc5[i][0] += nn_output[0];
512 solacc5[i][1] += nn_output[1];
514 test_size (j, psisol0REC.NoOfInput, "psisol0");
516 psisol0 (netprofin3, nn_output, 0);
517 solacc0[i][0] += nn_output[0];
518 solacc0[i][1] += nn_output[1];
522 for (i = 0; i < data->lens; i++) {
523 if (solacc25[i][0] > solacc25[i][1]) {
526 if (solacc25[i][1] > solacc25[i][0]) {
529 if (solacc5[i][0] > solacc5[i][1]) {
532 if (solacc5[i][1] > solacc5[i][0]) {
535 if (solacc0[i][0] > solacc0[i][1]) {
538 if (solacc0[i][1] > solacc0[i][0]) {
544 float seq2str[MAXSEQLEN][3];
547 /* This is for the PSIBLAST PSSM, for which there are two networks for each layer */
548 memset (psi2net, 0x00, MAXSEQLEN * 3);
550 for (t = 0; t < 2; t++) {
554 for (i = 0; i < data->lens; i++) {
558 doprofpsi (data, windows, i, psiar);
560 for (y = 0; y < windows; y++) {
562 for (z = 0; z < 20; z++)
563 netprofin[j++] = psiar[y][z];
568 test_size (j, pssma1REC.NoOfInput, "pssma1");
569 pssma1 (netprofin, nn_output, 0);
571 test_size (j, pssmb1REC.NoOfInput, "pssmb1");
572 pssmb1 (netprofin, nn_output, 0);
575 seq2str[i][0] = nn_output[0];
576 seq2str[i][1] = nn_output[1];
577 seq2str[i][2] = nn_output[2];
582 for (i = 0; i < data->lens; i++) {
585 float nn_input[19 * 3];
587 dowinsspsi (seq2str, windows, i, winarss, data);
589 for (y = 0; y < windows; y++) {
591 for (z = 0; z < 3; z++)
592 nn_input[j++] = winarss[y][z];
597 test_size (j, pssma2REC.NoOfInput, "pssma2");
598 pssma2 (nn_input, nn_output, 0);
600 test_size (j, pssmb2REC.NoOfInput, "pssmb2");
601 pssmb2 (nn_input, nn_output, 0);
604 psi2net[i][0] += nn_output[0];
605 psi2net[i][1] += nn_output[1];
606 psi2net[i][2] += nn_output[2];
609 for (i = 0; i < data->lens; i++) {
615 for (i = 0; i < data->lens; i++) {
620 #endif /* end of ifdef PSSM */
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.
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];
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;
642 if (finalout[i][0] > finalout[i][1] && finalout[i][0] > finalout[i][2])
644 if (finalout[i][1] > finalout[i][0] && finalout[i][1] > finalout[i][2])
646 if (finalout[i][2] > finalout[i][0] && finalout[i][2] > finalout[i][1])
649 int2pred (psi2fin, psi2let, data->lens);
650 #endif /* end PSSM */
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];
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;
666 if (finalout[i][0] > finalout[i][1] && finalout[i][0] > finalout[i][2])
668 if (finalout[i][1] > finalout[i][0] && finalout[i][1] > finalout[i][2])
670 if (finalout[i][2] > finalout[i][0] && finalout[i][2] > finalout[i][1])
673 int2pred (hmmfin, hmmlet, data->lens);
675 #endif /* end ifdef HMM */
677 /* When all data is available use the 'cons' network to decide all non-jury positions.
678 * Also calculate the confidence values */
681 for (j = 0; j < data->lens; j++) {
683 /* Is there disagreement? */
684 if (psi2let[j] == hmmlet[j]) { /* No - use the PSI networks */
686 consfin[j] = psi2fin[j];
688 confidence[j] = doconf (psi2net[j][0], psi2net[j][1], psi2net[j][2]);
690 } else { /* yes - take the mean propensities of the HMM and PSI networks */
694 finalout[j][0] = hmmnet[j][0];
695 finalout[j][1] = hmmnet[j][1];
696 finalout[j][2] = hmmnet[j][2];
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 */
703 confidence[j] = doconf (finalout[j][0], finalout[j][1], finalout[j][2]);
707 if (finalout[j][0] > finalout[j][1] && finalout[j][0] > finalout[j][2])
709 if (finalout[j][1] > finalout[j][0] && finalout[j][1] > finalout[j][2])
711 if (finalout[j][2] > finalout[j][0] && finalout[j][2] > finalout[j][1])
714 } /* agreement if/else block */
716 } /* end of sequence */
717 consfin[data->lens] = 25;
720 /* end of if (NOHMM==0 && NOPSI==0) */
721 /* Calculate the confidence values */
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];
727 confidence[i] = doconf ((hmmnet[i][0]), (hmmnet[i][1]), (hmmnet[i][2]));
729 consfin[data->lens] = 25;
733 if (1 == TESTOUTPUT) {
734 fprintf (stderr, "\n\nBoth PSIBLAST and HMM profiles were found\nAccuracy will average 81.6%%\n\n");
736 fprintf (stdout, "\n\nBoth PSIBLAST and HMM profiles were found\nAccuracy will average 81.6%%\n\n");
739 int2pred (consfin, finlet, data->lens);
740 finlet[data->lens] = '\0';
743 /* Remove unlikely secondary structure */
744 filter (finlet, letfilt);
745 filter (finlet, letfilt);
746 filter (finlet, letfilt);
748 memcpy (letfilt, finlet, sizeof (char) * MAXSEQLEN);
751 /* output predictions */
752 if (HUMAN) /* print human readable and quit */
753 print_human (data->lens, letfilt, confidence);
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]);
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]);
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]);
796 /*************************************************************************/
797 /* Calculate confidence as for PHD method (Rost and Sander, JMB, 1993, v232, p584) */
799 doconf (const float confa, const float confb, const float confc)
801 float whichout, maxout, maxnext, outconf;
802 whichout = outconf = maxnext = 0;
805 if (confa > confb && confa > confc) {
809 if (confb > confa && confb > confc) {
813 if (confc > confa && confc > confb) {
836 outconf = 10 * (maxout - maxnext);
843 /*************************************************************************/
845 StrReplStr (char *s1, char *s2, char *FromSStr, char *ToSStr)
852 while ((ChP2 = strstr (ChP1, FromSStr)) != NULL) {
854 Strncat (s1, ChP1, strlen (ChP1) - strlen (ChP2));
856 ChP1 = ChP2 + strlen (FromSStr);
862 /*************************************************************************/
864 Strncat (char *s1, char *s2, int len)
868 fprintf (stderr, "ERROR! Strncat error!");
872 if (s1 == NULL || s2 == NULL) {
873 fprintf (stderr, "ERROR! Strncat error!");
876 OrigLen = strlen (s1);
877 if (strncat (s1, s2, len) == NULL) {
878 fprintf (stderr, "ERROR! Strncat error!");
882 s1[OrigLen + len] = '\0';
887 /*************************************************************************/
888 /* Finds the data around the current position in a window */
890 doprofpsi (const alldata * data, const int win, const int curpos, float psiar[30][30])
894 for (i = (curpos - ((win - 1) / 2)); i <= (curpos + ((win - 1) / 2)); i++) {
896 for (k = 0; k < 20; k++) {
897 psiar[j][k] = data->psimat[i][k];
900 if (i < 0 || i >= data->lens)
901 memset (psiar[j], 0x00, sizeof (psiar[0][0]) * 30);
907 /*************************************************************************/
908 /* Converts an array of predicted structure to a character array, not a null terminated string */
910 int2pred (const int *pred, char *letpred, const int length)
914 for (i = 0; i < length; i++) {
925 /*case 25: letpred[i] = '\0'; break; */
927 fprintf (stderr, "ERROR! unknown secondary structure type '%i' at position %i in int2pred \n", pred[i], i);
934 /*************************************************************************/
936 doprofhmm (const alldata * data, const int win, const int curpos, float psiar[30][30])
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];
945 if (i < 0 || i >= data->lens)
946 for (k = 0; k < 24; k++)
952 /*************************************************************************/
953 /** print out a human readable version of the prediction and confidence **/
955 print_human (int length, char pred[MAXSEQLEN], float conf[MAXSEQLEN])
957 int i, j = 0, k; /* counters */
959 float *confout; /* array pointers to temp storage of output lines */
962 /* allocate memory to array pointers */
963 confout = malloc (BLOCK * sizeof (float));
964 if (confout == NULL) {
965 fprintf (stderr, "ERROR! Out of memory\n");
969 predout = malloc (BLOCK * sizeof (char));
970 if (predout == NULL) {
971 fprintf (stderr, "ERROR! Out of memory\n");
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 */
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]);
986 for (k = 0; k < BLOCK; ++k)
987 printf ("%1.0f", confout[k]);
989 j = 0; /* zero the block counter */
995 /* catch any trailing blocks that are smaller than required */
997 printf ("\n\npred: ");
998 for (k = 0; k < j; ++k)
999 printf ("%c", predout[k]);
1001 printf ("\nconf: ");
1002 for (k = 0; k < j; ++k)
1003 printf ("%1.0f", confout[k]);
1006 /* free memory and exit */
1009 exit (EXIT_SUCCESS);