--- /dev/null
+/************************************************************************
+ * Jnet - A consensus neural network protein *
+ * secondary structure prediction method *
+ * *
+ * Copyright 1999,2009 James Cuff, Jonathan Barber and Christian Cole *
+ * *
+ ************************************************************************
+
+
+-------------------------------------------------------------------------
+ This program is free software: you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ This program is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with this program. If not, see <http://www.gnu.org/licenses/>.
+-------------------------------------------------------------------------*/
+
+/* Fri Nov 15 15:28:09 GMT 2002 Bug fixes by JDB */
+/* WARNING: MUST BE PRESENTED WITH UPPERCASE DATA as lowercase 'n'
+ * has a special internal meaning in seq2int
+ * */
+
+ /*#define POST *//* Allow post processing of NN output */
+ /*#define FILTER *//* Allow string post processing of final output */
+
+/** The profile-specific defines below are not independent of each other.
+Especially check that the Jury set-up is correct! */
+#define HMM /* Do HMMer predictions */
+#define PSSM /* Do PSSM predicitons */
+ /*#define HMMONLY *//* For non-jury positions use the HMM prediction - only works when DOJURY is undefined */
+
+/* The allowed value of the integers representing residues */
+enum
+{
+ WINAR = 30, /* array size for winar */
+ PSILEN = 20, /* no. of elements per position in PSIBLAST PSSM */
+ PROLEN = 24, /* no. of elements per position for other profiles */
+ HELIX = 1, /* enums for sec. struct. states */
+ SHEET = 0,
+ COIL = 2,
+ MAXSEQNUM = 1000, /* Maximum number of sequences per alignment */
+ MAXSEQLEN = 5000, /* Maximum Sequence Length */
+ MAXBLOCK = 200 /* Maximum block size for human readable output */
+};
+
+#include <stdio.h>
+#include <string.h>
+#include <stdlib.h>
+#include <math.h>
+#include <ctype.h>
+#include <assert.h>
+
+#include "cmdline/cmdline.h"
+
+#include "pssma1.h" /* Uses information from data.psimat, (PSSM) - 100 hidden nodes */
+#include "pssma2.h" /* Second network for pssma1 - 100 hidden nodes */
+#include "pssmb1.h" /* Uses information from data.psimat, (PSSM) - 20 hidden nodes */
+#include "pssmb2.h" /* Second network for pssmb1 - 20 hidden nodes */
+
+#include "hmm1.h" /* First network for HMMER profile prediction */
+#include "hmm2.h" /* Second network for hmm1 */
+
+#include "hmmsol25.h"
+#include "hmmsol5.h"
+#include "hmmsol0.h"
+
+#include "psisol25.h"
+#include "psisol5.h"
+#include "psisol0.h"
+
+/* Main data structure for information about input data. Not all the fields are
+ * used */
+typedef struct alldata
+{
+ int lens; /* length of sequences */
+ int profmat[MAXSEQLEN][PROLEN]; /* BLOSUM frequencies from PSIBLAST alignment - unused but can't remove without segfault? */
+ float psimat[MAXSEQLEN][PSILEN]; /* PSIBLAST PSSM */
+ float hmmmat[MAXSEQLEN][PROLEN]; /* results from HMMER */
+}
+alldata;
+
+/* Global Vars - used to specify which input profiles have been provided */
+int NOPSI = 0;
+int NOHMM = 0;
+int TESTOUTPUT = 0;
+int HUMAN = 0; /* toggle human readable output */
+int BLOCK = 60; /* block size for human output */
+
+/* PSSM-specific functions */
+void readpsi (FILE * psifile, alldata * data);
+void doprofpsi (const alldata * data, const int win, const int curpos, float psiar[30][30]);
+
+/* HMMer-specific functions */
+void readhmm (FILE * hmmfile, alldata * data);
+void doprofhmm (const alldata * data, const int win, const int curpos, float psiar[30][30]);
+
+/* General or non-specific functions */
+void int2pred (const int *pred, char *letpred, const int length);
+float doconf (const float confa, const float confb, const float confc);
+void dowinsspsi (float seq2str[MAXSEQLEN][3], int win, int curpos, float winarss[30][4], const alldata * data);
+void dopred (const alldata * data);
+void test_size (int size, int required_size, char *name);
+void print_human (int length, char pred[MAXSEQLEN], float conf[MAXSEQLEN]);
+
+/* Prediction post-processing functions - not be used anymore, but retained for testing */
+/* Executing controlled by #define FILTER */
+void filter (char inp[MAXSEQLEN], char outp[MAXSEQLEN]);
+void StrReplStr (char *s1, char *s2, char *FromSStr, char *ToSStr);
+char *Strncat (char *s1, char *s2, int len);
+
+/*************************************************************/
+int
+main (int argc, char **argv)
+{
+ FILE *psifile, *hmmfile;
+ alldata *data;
+
+ /* Declared in cmdline/cmdline.h */
+ struct gengetopt_args_info args_info;
+
+ /* Set up the data structure for the input data */
+ if ((data = malloc (sizeof (alldata))) == NULL) {
+ fprintf (stderr, "ERROR! malloc failed at line %i\n", __LINE__);
+ exit (EXIT_FAILURE);
+ }
+
+ /* Get the command line arguments */
+ if (cmdline_parser (argc, argv, &args_info) != 0)
+ exit (EXIT_FAILURE);
+
+ /* check whether TEST STYLE OUTPUT has been requested */
+ if (args_info.test_flag) {
+ TESTOUTPUT = 1;
+ }
+
+ /* Print version info - CMDLINE_PARSER_PACKAGE CMDLINE_PARSER_VERSION declared in cmdline/cmdline.h */
+ if (1 == TESTOUTPUT) {
+ fprintf (stderr, "%s %s\n", CMDLINE_PARSER_PACKAGE, CMDLINE_PARSER_VERSION);
+ } else {
+ fprintf (stdout, "%s %s\n", CMDLINE_PARSER_PACKAGE, CMDLINE_PARSER_VERSION);
+ }
+
+ /* Read the HMMER file. Don't need to check if it is given as is a required option */
+ if ((hmmfile = fopen (args_info.hmmer_arg, "r")) == NULL) {
+ fprintf (stderr, "ERROR: Can't open HMMER file %s\n", args_info.hmmer_arg);
+ exit (EXIT_FAILURE);
+ } else {
+ if (1 == TESTOUTPUT) {
+ fprintf (stderr, "Found HMM profile data\n");
+ } else {
+ fprintf (stdout, "Found HMM profile data\n");
+ }
+ readhmm (hmmfile, data);
+ (void) fclose (hmmfile);
+ }
+
+ /* Read the PSIBLAST PSSM */
+ if (args_info.pssm_given) {
+ if ((psifile = fopen (args_info.pssm_arg, "r")) == NULL) {
+ fprintf (stderr, "ERROR: Can't open PSSM profile file %s\n", args_info.pssm_arg);
+ exit (EXIT_FAILURE);
+ } else {
+ if (1 == TESTOUTPUT) {
+ fprintf (stderr, "Found PSSM profile file\n");
+ } else {
+ fprintf (stdout, "Found PSSM profile file\n");
+ }
+ readpsi (psifile, data);
+ (void) fclose (psifile);
+ }
+ } else {
+ NOPSI = 1;
+ }
+
+ /* check whether human readable output is required (set by --human) */
+ if (args_info.human_given) {
+ HUMAN = 1;
+ if (args_info.human_orig != NULL) {
+ printf ("human: %d\n", args_info.human_arg);
+ printf ("human orig: %s\n", args_info.human_orig);
+ BLOCK = args_info.human_arg;
+ if (BLOCK > MAXBLOCK) {
+ fprintf (stderr, "ERROR - argument for --human is too large\n");
+ exit (EXIT_FAILURE);
+ }
+ }
+ }
+
+ /* Do the prediction */
+ if (1 == TESTOUTPUT) {
+ fprintf (stderr, "Running final predictions!\n");
+ } else {
+ fprintf (stdout, "Running final predictions!\n");
+ }
+ dopred (data);
+
+ free (data);
+
+ return EXIT_SUCCESS;
+}
+
+
+/*************************************************************************/
+/* Replaces particular runs of secondary structure with something else */
+void
+filter (char inp[MAXSEQLEN], char outp[MAXSEQLEN])
+{
+ int i;
+
+ struct filter_pair
+ {
+ char *orig;
+ char *final;
+ } pairs[] = {
+ {
+ "EHHHE", "EEEEE"}, {
+ "-HHH-", "HHHHH"}, {
+ "EHHH-", "EHHHH"}, {
+ "HHHE-", "HHH--"}, {
+ "-EHHH", "--HHH"}, {
+ "EHHE", "EEEE"}, {
+ "-HEH-", "-HHH-"}, {
+ "EHH-", "EEE-"}, {
+ "-HHE", "-EEE"}, {
+ "-HH-", "----"}, {
+ "HEEH", "EEEE"}, {
+ "-HE", "--E"}, {
+ "EH-", "E--"}, {
+ "-H-", "---"}, {
+ "HEH", "HHH"}, {
+ "-E-E-", "-EEE-"}, {
+ "E-E", "EEE"}, {
+ "H-H", "HHH"}, {
+ "EHE", "EEE"}, {
+ NULL, NULL}
+ };
+
+ for (i = 0; pairs[i].orig != NULL; i++) {
+ StrReplStr (outp, inp, pairs[i].orig, pairs[i].final);
+ strcpy (inp, outp);
+ }
+}
+
+/*************************************************************************/
+/* Reads in the PSIBLAST PSSM */
+void
+readpsi (FILE * psifile, alldata * data)
+{
+ int nlines = 0;
+
+ while ((getc (psifile)) != EOF) {
+ float tempval;
+ int i;
+ for (i = 0; i < 20; i++) {
+ int npar = fscanf (psifile, "%f", &tempval);
+ if (1 == npar) {
+ data->psimat[nlines][i] = tempval;
+ } else {
+ fprintf (stderr, "ERROR: Can't read a value from PSSM file");
+ data->psimat[nlines][i] = 0.;
+ }
+ }
+ nlines++;
+ }
+
+ /* check the data length matches that of the HMMer data */
+ if (data->lens != nlines - 1) {
+ fprintf (stderr, "ERROR! PSSM data length doesn't match that for HMMer. Are both input files from the same alignment/sequence?\n");
+ exit (EXIT_FAILURE);
+ }
+}
+
+
+/*************************************************************************/
+/* Reads in HMMER profile that was originally based upon ClustalW alignments */
+/* Additionally, use this data to define the length of the sequence in the absence of an alignment */
+void
+readhmm (FILE * hmmfile, alldata * data)
+{
+ int nlines = 0;
+
+ while ((getc (hmmfile)) != EOF) {
+ int i;
+
+ if (nlines > MAXSEQLEN) {
+ fprintf (stderr, "ERROR! the HMMER profile is too long. Can't cope... must quit.\n");
+ exit (EXIT_FAILURE);
+ }
+
+ for (i = 0; i < 24; i++) {
+ float tempval;
+ int npar = fscanf (hmmfile, "%f", &tempval);
+ if (1 == npar) {
+ data->hmmmat[nlines][i] = tempval;
+ } else {
+ fprintf (stderr, "ERROR: Can't read a value from HMMER profile file");
+ data->hmmmat[nlines][i] = 0.;
+ }
+ }
+ nlines++;
+ }
+ data->lens = nlines - 1; /* set the query sequence length */
+}
+
+/*************************************************************************/
+/* Function to test whether the data array for input
+ into the Neural Network function is of the expected size */
+void
+test_size (int size, int required_size, char *name)
+{
+ if (size != required_size) {
+ fprintf (stderr, "ERROR! input array not of correct size for %s (%i vs. %i)\n", name, size, required_size);
+ exit (EXIT_FAILURE);
+ }
+}
+
+
+
+/*************************************************************************/
+/***********************************************************************/
+/* Functions below here haven't been checked */
+void
+dowinsspsi (float seq2str[MAXSEQLEN][3], int win, int curpos, float winarss[30][4], const alldata * data)
+{
+ int i, j;
+
+ for (j = 0, i = (curpos - ((win - 1) / 2)); i <= (curpos + ((win - 1) / 2)); i++, j++) {
+
+ if (i >= 0 && i < data->lens) {
+ /* CC - changesd to make it more similar to network training */
+
+ /* for (k = 0; k < 3; k++)
+ winarss[j][k] = seq2str[i][k]; */
+
+ if ((seq2str[i][0] > seq2str[i][1]) && (seq2str[i][0] > seq2str[i][2])) {
+ winarss[j][0] = 1.0;
+ winarss[j][1] = 0.0;
+ winarss[j][2] = 0.0;
+ } else if ((seq2str[i][1] > seq2str[i][0]) && (seq2str[i][1] > seq2str[i][2])) {
+ winarss[j][0] = 0.0;
+ winarss[j][1] = 1.0;
+ winarss[j][2] = 0.0;
+ } else if ((seq2str[i][2] > seq2str[i][0]) && (seq2str[i][2] > seq2str[i][1])) {
+ winarss[j][0] = 0.0;
+ winarss[j][1] = 0.0;
+ winarss[j][2] = 1.0;
+ } else {
+ winarss[j][0] = 0.0;
+ winarss[j][1] = 0.0;
+ winarss[j][2] = 1.0;
+ }
+ } else if (i < 0 || i >= data->lens) {
+ /* Setting winarss[j][2] to 1 increases q3 accuracy */
+ /* CC - does it really??? */
+ winarss[j][0] = 0.0;
+ winarss[j][1] = 0.0;
+ winarss[j][2] = 0.0;
+ }
+ }
+}
+
+/*************************************************************************/
+/* Structure to hold data from networks */
+typedef struct netdata
+{
+ int ps1;
+ float *netin;
+ float *netprofin;
+} netdata;
+
+/*************************************************************************/
+/* main prediction routine */
+void
+dopred (const alldata * data)
+{
+ extern int NOHMM, NOPSI;
+
+ float psi2net[MAXSEQLEN][3];
+ float hmmnet[MAXSEQLEN][3];
+
+ float finalout[MAXSEQLEN][3];
+ int psi2fin[MAXSEQLEN];
+ int hmmfin[MAXSEQLEN];
+ int consfin[MAXSEQLEN];
+
+ char psi2let[MAXSEQLEN];
+ char hmmlet[MAXSEQLEN];
+ char finlet[MAXSEQLEN];
+
+ char sollet25[MAXSEQLEN];
+ char sollet5[MAXSEQLEN];
+ char sollet0[MAXSEQLEN];
+
+ float netprofin3[500];
+ float confidence[MAXSEQLEN]; /* Array for holding the confidence value calculated by doconf() */
+ float psiar[WINAR][WINAR];
+ int winar2[WINAR][WINAR];
+
+ float solacc25[MAXSEQLEN][2];
+ float solacc5[MAXSEQLEN][2];
+ float solacc0[MAXSEQLEN][2];
+
+ float winarss[WINAR][4];
+
+ char letfilt[MAXSEQLEN];
+
+ int i, t;
+ char jury[MAXSEQLEN];
+
+ float netprofin[MAXSEQLEN];
+ memset (netprofin, 0x00, sizeof (float) * MAXSEQLEN);
+ memset (winar2, 0x00, sizeof (winar2[0][0]) * WINAR * WINAR);
+
+ letfilt[0] = '\0';
+
+#ifdef HMM
+ {
+ float seq2str[MAXSEQLEN][3];
+ int i, windows = 17;
+ for (i = 0; i < data->lens; i++) {
+ float nn_output[3];
+
+ int y, j = 0;
+ doprofhmm (data, windows, i, psiar);
+
+ for (y = 0; y < windows; y++) {
+ int z;
+ for (z = 0; z < PROLEN; z++)
+ netprofin3[j++] = psiar[y][z];
+ }
+
+ test_size (j, hmm1REC.NoOfInput, "hmm1");
+
+ hmm1 (netprofin3, nn_output, 0);
+ seq2str[i][0] = nn_output[0];
+ seq2str[i][1] = nn_output[1];
+ seq2str[i][2] = nn_output[2];
+ }
+
+ windows = 19;
+ for (i = 0; i < data->lens; i++) {
+ int y, j = 0;
+ float nn_output[3];
+ float nn_input[19 * 3];
+
+ dowinsspsi (seq2str, windows, i, winarss, data);
+
+ for (y = 0; y < windows; y++) {
+ int z;
+ for (z = 0; z < 3; z++)
+ nn_input[j++] = winarss[y][z];
+ }
+
+ test_size (j, hmm2REC.NoOfInput, "hmm2");
+
+ hmm2 (nn_input, nn_output, 0);
+ hmmnet[i][0] = nn_output[0];
+ hmmnet[i][1] = nn_output[1];
+ hmmnet[i][2] = nn_output[2];
+ }
+ }
+#endif /* end ifdef HMM */
+
+ if (NOHMM == 0) {
+ int i, windows = 17;
+ for (i = 0; i < data->lens; i++) {
+ float nn_output[3];
+ int y, j = 0;
+ doprofhmm (data, windows, i, psiar);
+
+ for (y = 0; y < windows; y++) {
+ int z;
+ for (z = 0; z < 24; z++)
+ netprofin3[j++] = psiar[y][z];
+ }
+
+ test_size (j, hmmsol25REC.NoOfInput, "hmmsol25");
+
+ hmmsol25 (netprofin3, nn_output, 0);
+ solacc25[i][0] = nn_output[0];
+ solacc25[i][1] = nn_output[1];
+
+ test_size (j, hmmsol5REC.NoOfInput, "hmmsol5");
+
+ hmmsol5 (netprofin3, nn_output, 0);
+ solacc5[i][0] = nn_output[0];
+ solacc5[i][1] = nn_output[1];
+
+ test_size (j, hmmsol0REC.NoOfInput, "hmmsol0");
+
+ hmmsol0 (netprofin3, nn_output, 0);
+ solacc0[i][0] = nn_output[0];
+ solacc0[i][1] = nn_output[1];
+ }
+ }
+
+ if (NOPSI == 0) {
+ int i, windows = 17;
+ for (i = 0; i < data->lens; i++) {
+ float nn_output[3];
+ int y, j = 0;
+ doprofpsi (data, windows, i, psiar);
+
+ for (y = 0; y < windows; y++) {
+ int z;
+ for (z = 0; z < 20; z++)
+ netprofin3[j++] = psiar[y][z];
+ }
+
+ test_size (j, psisol25REC.NoOfInput, "psisol25");
+
+ psisol25 (netprofin3, nn_output, 0);
+ solacc25[i][0] += nn_output[0];
+ solacc25[i][1] += nn_output[1];
+
+ test_size (j, psisol5REC.NoOfInput, "psisol5");
+
+ psisol5 (netprofin3, nn_output, 0);
+ solacc5[i][0] += nn_output[0];
+ solacc5[i][1] += nn_output[1];
+
+ test_size (j, psisol0REC.NoOfInput, "psisol0");
+
+ psisol0 (netprofin3, nn_output, 0);
+ solacc0[i][0] += nn_output[0];
+ solacc0[i][1] += nn_output[1];
+ }
+ }
+
+ for (i = 0; i < data->lens; i++) {
+ if (solacc25[i][0] > solacc25[i][1]) {
+ sollet25[i] = '-';
+ }
+ if (solacc25[i][1] > solacc25[i][0]) {
+ sollet25[i] = 'B';
+ }
+ if (solacc5[i][0] > solacc5[i][1]) {
+ sollet5[i] = '-';
+ }
+ if (solacc5[i][1] > solacc5[i][0]) {
+ sollet5[i] = 'B';
+ }
+ if (solacc0[i][0] > solacc0[i][1]) {
+ sollet0[i] = '-';
+ }
+ if (solacc0[i][1] > solacc0[i][0]) {
+ sollet0[i] = 'B';
+ }
+ }
+
+ if (NOPSI == 0) {
+ float seq2str[MAXSEQLEN][3];
+ int i, windows = 17;
+
+ /* This is for the PSIBLAST PSSM, for which there are two networks for each layer */
+ memset (psi2net, 0x00, MAXSEQLEN * 3);
+#ifdef PSSM
+ for (t = 0; t < 2; t++) {
+
+ /* First network */
+ windows = 17;
+ for (i = 0; i < data->lens; i++) {
+ int y, j = 0;
+ float nn_output[3];
+
+ doprofpsi (data, windows, i, psiar);
+
+ for (y = 0; y < windows; y++) {
+ int z;
+ for (z = 0; z < 20; z++)
+ netprofin[j++] = psiar[y][z];
+ }
+
+
+ if (t == 0) {
+ test_size (j, pssma1REC.NoOfInput, "pssma1");
+ pssma1 (netprofin, nn_output, 0);
+ } else if (t == 1) {
+ test_size (j, pssmb1REC.NoOfInput, "pssmb1");
+ pssmb1 (netprofin, nn_output, 0);
+ }
+
+ seq2str[i][0] = nn_output[0];
+ seq2str[i][1] = nn_output[1];
+ seq2str[i][2] = nn_output[2];
+ }
+
+ /* Second network */
+ windows = 19;
+ for (i = 0; i < data->lens; i++) {
+ int y, j = 0;
+ float nn_output[3];
+ float nn_input[19 * 3];
+
+ dowinsspsi (seq2str, windows, i, winarss, data);
+
+ for (y = 0; y < windows; y++) {
+ int z;
+ for (z = 0; z < 3; z++)
+ nn_input[j++] = winarss[y][z];
+ }
+
+
+ if (t == 0) {
+ test_size (j, pssma2REC.NoOfInput, "pssma2");
+ pssma2 (nn_input, nn_output, 0);
+ } else if (t == 1) {
+ test_size (j, pssmb2REC.NoOfInput, "pssmb2");
+ pssmb2 (nn_input, nn_output, 0);
+ }
+
+ psi2net[i][0] += nn_output[0];
+ psi2net[i][1] += nn_output[1];
+ psi2net[i][2] += nn_output[2];
+ }
+ }
+ for (i = 0; i < data->lens; i++) {
+ psi2net[i][0] /= 2;
+ psi2net[i][1] /= 2;
+ psi2net[i][2] /= 2;
+ }
+#else
+ for (i = 0; i < data->lens; i++) {
+ psi2net[i][0] = 0;
+ psi2net[i][1] = 0;
+ psi2net[i][2] = 0;
+ }
+#endif /* end of ifdef PSSM */
+ }
+
+ /* Work out the final output values for the various outputs
+ * There is some smoothing going on at the ends of the sequences
+ * Also, some of these are being done even if there is no data in them.
+ * */
+
+ if (NOPSI == 0) {
+#ifdef PSSM
+ for (i = 0; i < data->lens; i++) {
+ finalout[i][0] = psi2net[i][0];
+ finalout[i][1] = psi2net[i][1];
+ finalout[i][2] = psi2net[i][2];
+
+#ifdef POST
+ if (i <= 4)
+ finalout[i][2] += (5 - i) * 0.2;
+ else if (i >= (data->lens - 5))
+ finalout[i][2] += (5 - ((data->lens - 1) - i)) * 0.2;
+#endif
+
+ if (finalout[i][0] > finalout[i][1] && finalout[i][0] > finalout[i][2])
+ psi2fin[i] = SHEET;
+ if (finalout[i][1] > finalout[i][0] && finalout[i][1] > finalout[i][2])
+ psi2fin[i] = HELIX;
+ if (finalout[i][2] > finalout[i][0] && finalout[i][2] > finalout[i][1])
+ psi2fin[i] = COIL;
+ }
+ int2pred (psi2fin, psi2let, data->lens);
+#endif /* end PSSM */
+ }
+#ifdef HMM
+ if (NOHMM == 0) {
+ for (i = 0; i < data->lens; i++) {
+ finalout[i][0] = hmmnet[i][0];
+ finalout[i][1] = hmmnet[i][1];
+ finalout[i][2] = hmmnet[i][2];
+
+#ifdef POST
+ if (i <= 4)
+ finalout[i][2] += (5 - i) * 0.2;
+ else if (i >= (data->lens - 5))
+ finalout[i][2] += (5 - ((data->lens - 1) - i)) * 0.2;
+#endif
+
+ if (finalout[i][0] > finalout[i][1] && finalout[i][0] > finalout[i][2])
+ hmmfin[i] = SHEET;
+ if (finalout[i][1] > finalout[i][0] && finalout[i][1] > finalout[i][2])
+ hmmfin[i] = HELIX;
+ if (finalout[i][2] > finalout[i][0] && finalout[i][2] > finalout[i][1])
+ hmmfin[i] = COIL;
+ }
+ int2pred (hmmfin, hmmlet, data->lens);
+ }
+#endif /* end ifdef HMM */
+
+ /* When all data is available use the 'cons' network to decide all non-jury positions.
+ * Also calculate the confidence values */
+ if (NOPSI == 0) {
+ int j;
+ for (j = 0; j < data->lens; j++) {
+
+ /* Is there disagreement? */
+ if (psi2let[j] == hmmlet[j]) { /* No - use the PSI networks */
+ jury[j] = ' ';
+ consfin[j] = psi2fin[j];
+
+ confidence[j] = doconf (psi2net[j][0], psi2net[j][1], psi2net[j][2]);
+
+ } else { /* yes - take the mean propensities of the HMM and PSI networks */
+ jury[j] = '*';
+
+#ifdef HMMONLY
+ finalout[j][0] = hmmnet[j][0];
+ finalout[j][1] = hmmnet[j][1];
+ finalout[j][2] = hmmnet[j][2];
+#else
+ finalout[j][0] = (hmmnet[j][0] + psi2net[j][0]) / 2;
+ finalout[j][1] = (hmmnet[j][1] + psi2net[j][1]) / 2;
+ finalout[j][2] = (hmmnet[j][2] + psi2net[j][2]) / 2;
+#endif /* end HMMONLY */
+
+ confidence[j] = doconf (finalout[j][0], finalout[j][1], finalout[j][2]);
+
+ consfin[j] = COIL;
+
+ if (finalout[j][0] > finalout[j][1] && finalout[j][0] > finalout[j][2])
+ consfin[j] = SHEET;
+ if (finalout[j][1] > finalout[j][0] && finalout[j][1] > finalout[j][2])
+ consfin[j] = HELIX;
+ if (finalout[j][2] > finalout[j][0] && finalout[j][2] > finalout[j][1])
+ consfin[j] = COIL;
+
+ } /* agreement if/else block */
+
+ } /* end of sequence */
+ consfin[data->lens] = 25;
+ }
+
+ /* end of if (NOHMM==0 && NOPSI==0) */
+ /* Calculate the confidence values */
+ if (NOPSI == 1) {
+ fprintf (stdout, "\n\nWARNING!: Only using the HMM profile\nAccuracy will average 79.6%%\n\n");
+ for (i = 0; i < data->lens; i++) {
+ consfin[i] = hmmfin[i];
+
+ confidence[i] = doconf ((hmmnet[i][0]), (hmmnet[i][1]), (hmmnet[i][2]));
+ }
+ consfin[data->lens] = 25;
+ }
+
+ if (NOPSI == 0) {
+ if (1 == TESTOUTPUT) {
+ fprintf (stderr, "\n\nBoth PSIBLAST and HMM profiles were found\nAccuracy will average 81.6%%\n\n");
+ } else {
+ fprintf (stdout, "\n\nBoth PSIBLAST and HMM profiles were found\nAccuracy will average 81.6%%\n\n");
+ }
+ }
+ int2pred (consfin, finlet, data->lens);
+ finlet[data->lens] = '\0';
+
+#ifdef FILTER
+ /* Remove unlikely secondary structure */
+ filter (finlet, letfilt);
+ filter (finlet, letfilt);
+ filter (finlet, letfilt);
+#else
+ memcpy (letfilt, finlet, sizeof (char) * MAXSEQLEN);
+#endif
+
+ /* output predictions */
+ if (HUMAN) /* print human readable and quit */
+ print_human (data->lens, letfilt, confidence);
+
+ printf ("\njnetpred:");
+ for (i = 0; i < data->lens; i++)
+ printf ("%c,", letfilt[i]);
+ printf ("\nJNETCONF:");
+ for (i = 0; i < data->lens; i++)
+ printf ("%1.0f,", confidence[i]);
+ printf ("\nJNETSOL25:");
+ for (i = 0; i < data->lens; i++)
+ printf ("%c,", sollet25[i]);
+ printf ("\nJNETSOL5:");
+ for (i = 0; i < data->lens; i++)
+ printf ("%c,", sollet5[i]);
+ printf ("\nJNETSOL0:");
+ for (i = 0; i < data->lens; i++)
+ printf ("%c,", sollet0[i]);
+ printf ("\nJNETHMM:");
+ for (i = 0; i < data->lens; i++)
+ printf ("%c,", hmmlet[i]);
+
+ if (NOPSI == 0) {
+ printf ("\nJNETPSSM:");
+ for (i = 0; i < data->lens; i++)
+ printf ("%c,", psi2let[i]);
+ printf ("\nJNETJURY:");
+ for (i = 0; i < data->lens; i++)
+ printf ("%c,", jury[i]);
+ }
+
+ printf ("\nJNETPROPE:");
+ for (i = 0; i < data->lens; i++)
+ printf ("%.4f,", hmmnet[i][0]);
+ printf ("\nJNETPROPH:");
+ for (i = 0; i < data->lens; i++)
+ printf ("%.4f,", hmmnet[i][1]);
+ printf ("\nJNETPROPC:");
+ for (i = 0; i < data->lens; i++)
+ printf ("%.4f,", hmmnet[i][2]);
+ printf ("\n");
+}
+
+
+/*************************************************************************/
+/* Calculate confidence as for PHD method (Rost and Sander, JMB, 1993, v232, p584) */
+float
+doconf (const float confa, const float confb, const float confc)
+{
+ float whichout, maxout, maxnext, outconf;
+ whichout = outconf = maxnext = 0;
+ maxout = confc;
+
+ if (confa > confb && confa > confc) {
+ whichout = 0;
+ maxout = confa;
+ }
+ if (confb > confa && confb > confc) {
+ whichout = 1;
+ maxout = confb;
+ }
+ if (confc > confa && confc > confb) {
+ whichout = 2;
+ maxout = confc;
+ }
+
+ if (whichout == 0) {
+ if (confb > confc)
+ maxnext = confb;
+ if (confc > confb)
+ maxnext = confc;
+ }
+ if (whichout == 1) {
+ if (confc > confa)
+ maxnext = confc;
+ if (confa > confc)
+ maxnext = confa;
+ }
+ if (whichout == 2) {
+ if (confb > confa)
+ maxnext = confb;
+ if (confa > confb)
+ maxnext = confa;
+ }
+ outconf = 10 * (maxout - maxnext);
+ if (outconf > 9)
+ outconf = 9;
+
+ return outconf;
+}
+
+/*************************************************************************/
+void
+StrReplStr (char *s1, char *s2, char *FromSStr, char *ToSStr)
+{
+ char *ChP1, *ChP2;
+
+ s1[0] = '\0';
+ ChP1 = s2;
+
+ while ((ChP2 = strstr (ChP1, FromSStr)) != NULL) {
+ if (ChP1 != ChP2)
+ Strncat (s1, ChP1, strlen (ChP1) - strlen (ChP2));
+ strcat (s1, ToSStr);
+ ChP1 = ChP2 + strlen (FromSStr);
+ }
+ strcat (s1, ChP1);
+ return;
+}
+
+/*************************************************************************/
+char *
+Strncat (char *s1, char *s2, int len)
+{
+ int OrigLen = 0;
+ if (len == 0) {
+ fprintf (stderr, "ERROR! Strncat error!");
+ return s1;
+ }
+
+ if (s1 == NULL || s2 == NULL) {
+ fprintf (stderr, "ERROR! Strncat error!");
+ return NULL;
+ }
+ OrigLen = strlen (s1);
+ if (strncat (s1, s2, len) == NULL) {
+ fprintf (stderr, "ERROR! Strncat error!");
+ return NULL;
+ }
+
+ s1[OrigLen + len] = '\0';
+
+ return s1;
+}
+
+/*************************************************************************/
+/* Finds the data around the current position in a window */
+void
+doprofpsi (const alldata * data, const int win, const int curpos, float psiar[30][30])
+{
+ int i, j = 0;
+
+ for (i = (curpos - ((win - 1) / 2)); i <= (curpos + ((win - 1) / 2)); i++) {
+ int k;
+ for (k = 0; k < 20; k++) {
+ psiar[j][k] = data->psimat[i][k];
+ }
+
+ if (i < 0 || i >= data->lens)
+ memset (psiar[j], 0x00, sizeof (psiar[0][0]) * 30);
+
+ j++;
+ }
+}
+
+/*************************************************************************/
+/* Converts an array of predicted structure to a character array, not a null terminated string */
+void
+int2pred (const int *pred, char *letpred, const int length)
+{
+ int i;
+
+ for (i = 0; i < length; i++) {
+ switch (pred[i]) {
+ case HELIX:
+ letpred[i] = 'H';
+ break;
+ case SHEET:
+ letpred[i] = 'E';
+ break;
+ case COIL:
+ letpred[i] = '-';
+ break;
+ /*case 25: letpred[i] = '\0'; break; */
+ default:
+ fprintf (stderr, "ERROR! unknown secondary structure type '%i' at position %i in int2pred \n", pred[i], i);
+ exit (EXIT_FAILURE);
+ break;
+ }
+ }
+}
+
+/*************************************************************************/
+void
+doprofhmm (const alldata * data, const int win, const int curpos, float psiar[30][30])
+{
+ int i, j, k;
+ j = k = i = 0;
+
+ for (i = (curpos - ((win - 1) / 2)); i <= (curpos + ((win - 1) / 2)); i++) {
+ for (k = 0; k < 24; k++)
+ psiar[j][k] = data->hmmmat[i][k];
+
+ if (i < 0 || i >= data->lens)
+ for (k = 0; k < 24; k++)
+ psiar[j][k] = 0;
+ j++;
+ }
+}
+
+/*************************************************************************/
+/* print out a human readable version of the prediction and confidence */
+void
+print_human (int length, char pred[MAXSEQLEN], float conf[MAXSEQLEN])
+{
+ int i, j = 0, k;
+ float *confout; /* array pointers to temp storage of output lines */
+ char *predout;
+
+ /* allocate memory to array pointers */
+ confout = malloc (BLOCK * sizeof (float));
+ if (confout == NULL) {
+ fprintf (stderr, "ERROR! Out of memory\n");
+ exit (EXIT_FAILURE);
+ }
+
+ predout = malloc (BLOCK * sizeof (char));
+ if (predout == NULL) {
+ fprintf (stderr, "ERROR! Out of memory\n");
+ exit (EXIT_FAILURE);
+ }
+
+ for (i = 0; i < length; i++) {
+ confout[j] = conf[i]; /* store current conf value */
+ predout[j] = pred[i]; /* store current prediction value */
+
+ /* when reached end of block, print stuff out */
+ if (i > 1 && i % BLOCK == 0) {
+ printf ("\n\npred: ");
+ for (k = 0; k < BLOCK; ++k)
+ printf ("%c", predout[k]);
+
+ printf ("\nconf: ");
+ for (k = 0; k < BLOCK; ++k)
+ printf ("%1.0f", confout[k]);
+
+ j = 0; /* zero the block counter */
+ }
+
+ ++j;
+ }
+
+ /* catch any trailing blocks that are smaller than required */
+ if (j)
+ printf ("\n\npred: ");
+ for (k = 0; k < j; ++k)
+ printf ("%c", predout[k]);
+
+ printf ("\nconf: ");
+ for (k = 0; k < j; ++k)
+ printf ("%1.0f", confout[k]);
+ printf ("\n");
+
+ /* free memory and exit */
+ free (confout);
+ free (predout);
+ exit (EXIT_SUCCESS);
+}
git://source.jalview.org / jpred.git / blobdiff