--- /dev/null
+/* Interface to sc2seq created by David Wilson and modified by Ethan Wolf */
+/* Other interfaces by Ethan Wolf 1995. */
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include <string.h>
+#include "scio.h"
+#include "switches.h"
+#include "likelihood.h"
+#include "sc.h"
+#include "scconst.h"
+#include "options.h"
+#include "compute_like.h"
+#include "interface.h"
+#include "stats.h"
+
+#define PIR_LIKES "/tmp/pir.likes"
+#define PS_OF_LIKELIHOOD "~/likelihood.ps"
+#define PAIR_PIR_TXT "/tmp/pir1.txt"
+#define SINGLE_PIR_TXT "/tmp/pir2.txt"
+
+double many_pprobs[MAX_TABLE_NUMBER][NUM_RES_TYPE][POSNUM],
+ many_pprobp[MAX_TABLE_NUMBER][NUM_RES_TYPE][NUM_RES_TYPE][POSNUM][POSNUM];
+
+
+double gprobs[NUM_RES_TYPE], gprobp[NUM_RES_TYPE][NUM_RES_TYPE][POSNUM];
+
+double (*pprobs)[POSNUM];
+double (*pprobp)[NUM_RES_TYPE][POSNUM][POSNUM];
+
+
+long pfreqs[MAX_TABLE_NUMBER][AANUM][POSNUM];
+long ptotals[MAX_TABLE_NUMBER][POSNUM];
+long pfreqp[MAX_TABLE_NUMBER][AANUM][AANUM][POSNUM][POSNUM];
+long ptotalp[MAX_TABLE_NUMBER][POSNUM][POSNUM];
+
+
+#define MALLOC(NUMBER,TYPE) ((TYPE*)Malloc((NUMBER)*sizeof(TYPE)))
+
+
+
+void *Malloc (n)
+ int n;
+{
+ void *ptr;
+ ptr = (void*)malloc(n);
+ if (ptr) return ptr;
+ fprintf(stderr,"Memory shortage!!\n");
+ exit(-1);
+}
+
+
+
+
+void switch_tables(int table)
+{
+ extern int WINDOW, window_length[MAX_TABLE_NUMBER];
+
+ pprobs= many_pprobs[table];
+ pprobp= many_pprobp[table];
+ WINDOW = window_length[table];
+}
+
+
+
+/******** Used to subtract off a sequence from the table. ***/
+
+void recalc_prob(Sequence sequence, int input_table, int mode)
+{
+ int ProlineFreeWindow = mode & PROLINE_FREE_MODE;
+
+ /* Remove the sequence from the positive table before computing on it. */
+ add_seq2(pfreqs[input_table], ptotals[input_table], pfreqp[input_table],
+ ptotalp[input_table],
+ sequence.reg, sequence.seq, sequence.seqlen, -1);
+
+ if (!(mode & WEIGHTED_PROBS))
+ estimate_database_probs2(input_table, pfreqs[input_table],
+ ptotals[input_table], many_pprobs,
+ pfreqp[input_table], ptotalp[input_table],
+ many_pprobp, ProlineFreeWindow);
+/** Add seq back into the frequence counts. **/
+ add_seq2(pfreqs[input_table], ptotals[input_table], pfreqp[input_table],
+ ptotalp[input_table],
+ sequence.reg, sequence.seq, sequence.seqlen, 1);
+
+ #ifdef AVG_WEIRD_PROBS
+ calc_weird_pprob_avg(many_pprobp[input_table], many_pprobs[input_table]);
+ #else
+ calc_weird_pprob_sum(many_pprobp[input_table], many_pprobs[input_table]);
+ #endif
+
+}
+
+
+
+
+
+
+
+void PairCoilData (FILE *fgin, FILE **fpin,
+ int ProlineFree, int table_to_remove_from,
+ char *input_filename, int mode)
+{
+ long gfreqs[AANUM]; long gtotals;
+ long gfreqp[AANUM][AANUM][POSNUM]; long gtotalp[POSNUM];
+
+ extern long pfreqs[MAX_TABLE_NUMBER][AANUM][POSNUM];
+ extern long ptotals[MAX_TABLE_NUMBER][POSNUM];
+ extern long pfreqp[MAX_TABLE_NUMBER][AANUM][AANUM][POSNUM][POSNUM];
+ extern long ptotalp[MAX_TABLE_NUMBER][POSNUM][POSNUM];
+ extern double SCALE0;
+
+
+ int i=0;
+
+ if (fgin) {
+ readgenfreq(fgin, >otals, gfreqs, gtotalp, gfreqp, stderr);
+ calcgprobs(gfreqs, gtotals, gprobs);
+ calcgprobp(gfreqp, gtotalp, gprobp);
+ }
+ else setgprob_to_uniform(gprobs, gprobp);
+
+ #ifdef AVG_WEIRD_PROBS
+ calc_weird_gprob_avg(gprobp, gprobs);
+ #else
+ calc_weird_gprob_sum(gprobp, gprobs);
+ #endif
+
+
+ for (i=0; i<NUMBER_TABLES; i++) { /* Go through all input tables. */
+ if (fpin[i] != NULL) { /* Calculate probabilities for that table. */
+ readposfreq(fpin[i],ptotals[i],pfreqs[i], ptotalp[i], pfreqp[i], stderr);
+
+ if (i == table_to_remove_from)
+ remove_all_input_seq_from_table(pfreqs[table_to_remove_from],
+ ptotals[table_to_remove_from],
+ pfreqp[table_to_remove_from],
+ ptotalp[table_to_remove_from], input_filename, -1);
+
+ estimate_database_probs2(i, pfreqs[i], ptotals[i], many_pprobs,
+ pfreqp[i], ptotalp[i], many_pprobp,
+ ProlineFree);
+ }
+ else {
+ fprintf(stderr,"\nCould not make input table number %d.",i);
+ exit(-1);
+ }
+
+
+#ifdef AVG_WEIRD_PROBS /* make table. */
+ calc_weird_pprob_avg(many_pprobp[i], many_pprobs[i]);
+#else
+ calc_weird_pprob_sum(many_pprobp[i], many_pprobs[i]);
+#endif
+
+ }
+
+ switch_tables(0);
+
+}
+
+
+NEWCOILSinit ()
+{
+ initialize_relative();
+}
+
+#define NEGCCRATIO 30
+#define NEGMEAN 0.77
+#define NEGSD 0.20
+#define CCMEAN 1.63
+#define CCSD 0.24
+#define SQRT2PI 2.5066282746310002
+
+
+
+/** To compute the paircoil likelihood line. */
+
+
+void likelihood_from_txt(char *txt_filename, double *m, double *b1)
+{
+ FILE *txt_file;
+ int mode2=0;
+ double percent_coil=PERCENT_COILED;
+ double score1;
+ int filenum=1;
+ long n[2]={0,0}, n_halfgauss[2]={0,0};
+ double sum[2]={0,0}, sum2[2]={0,0};
+ double a[2],b[2]; /* ax+b is the likelihood line approximation */
+ double mean[2], sd[2], down[2];
+ long _counts[2][-4*MINSCORE +1];
+ double _means[2][-4*MINSCORE+1];
+ long *counts[2];
+ double *means[2];
+ int i, maxplaces[2]= {2*MINSCORE, 2*MINSCORE}, minplaces[2]={-2*MINSCORE,-2*MINSCORE};
+ int minplace, maxplace;
+ int number_strange=0, number_high=0;
+ FILE *devnull, *outfile, *outpsfile;
+ int place;
+
+ counts[0] = &(_counts[0][-2*MINSCORE]);
+ counts[1] = &(_counts[1][-2*MINSCORE]);
+ means[0] = &(_means[0][-2*MINSCORE]);
+ means[1] = &(_means[1][-2*MINSCORE]);
+
+
+ mode2 |= RATIO_TRAP_METHOD;
+ mode2 |= PERCENT_COIL_MODE;
+
+
+
+ for (i= 2*MINSCORE; i<=-2*MINSCORE; ++i) {
+ counts[0][i] = 0;
+ means[0][i]=0;
+ }
+
+
+ txt_file = sopen (txt_filename, "r");
+ while (fread(&score1,sizeof(double),1,txt_file)) {
+ if (score1<MINSCORE || score1>-(MINSCORE)) {
+/*
+ fprintf(stderr,"Strange score: %lf\n",score1);
+ exit(-1);
+*/
+ number_strange++;
+ }
+ place = ( (int)((score1-(MINSCORE))*2 + .5)) + 2*(MINSCORE);
+ if (place > -2*MINSCORE) {
+ fprintf(stderr, "Out of bounds: %d.\n", place);
+ fprintf(stderr, "Score: %lf\n",score1);
+ /* exit(-1); */
+ counts[0][-2*MINSCORE]++;
+ means[0][-2*MINSCORE]++;
+ number_high++;
+ maxplaces[0]=-2*MINSCORE;
+ }
+ else if (place < 2*MINSCORE) {
+ fprintf(stderr, "Out of bounds: %d.\n", place);
+ fprintf(stderr, "Score: %lf\n",score1);
+ counts[0][2*MINSCORE]++;
+ means[0][2*MINSCORE]++;
+ minplaces[0]=2*MINSCORE;
+ /* exit(-1); */
+ }
+
+ else {
+ counts[0][place]++;
+ means[0][place] += score1;
+ if (place>maxplaces[0]) maxplaces[0]=place;
+ if (place<minplaces[0]) minplaces[0]=place;
+ ++n[0];
+ }
+ }
+
+ fprintf(stderr,"\n%d out of bounds scores > %d were scored as %d\n",number_high,-(MINSCORE),-(MINSCORE));
+ fprintf(stderr,"\n%d strange scores < %d were scored as %d\n",number_strange-number_high,MINSCORE, MINSCORE);
+ sclose(txt_file);
+
+ /** Compute the means of the scores in each box. **/
+ for (i= minplaces[0]; i<=maxplaces[0]; ++i)
+ if (counts[0][i] != 0)
+ means[0][i] /= counts[0][i]; /** Each bucket gets average score*/
+
+
+ /*** We need a mean***/
+
+ mean[0]= percent_mean(counts[0], means[0], percent_coil,minplaces[0],
+ maxplaces[0],n[0]);
+
+
+ /*** Calculate gaussian fit to scores below the mean. **/
+ txt_file=sopen(txt_filename, "r");
+ while (fread(&score1,sizeof(double),1,txt_file))
+ if ( (score1 > MINSCORE) && (score1 < mean[0])) {
+ ++n_halfgauss[0];
+ sum2[0] += (score1-mean[0])*(score1-mean[0]);
+ }
+ sd[0] =sqrt(sum2[0]/n_halfgauss[0]);
+ down[0]= (double) 2*n_halfgauss[0]/n[0];
+
+ fprintf(stderr,"\n mean=%lf, sd=%lf, down= %lf", mean[0],sd[0],down[0]);
+ sclose(txt_file);
+
+ minplace=minplaces[0]; maxplace=maxplaces[0];
+
+ devnull=sopen("/dev/null","r+");
+ outfile=sopen(PIR_LIKES,"w");
+ outpsfile=sopen(PS_OF_LIKELIHOOD,"w");
+/* outfile=sopen("~/pir.out","w"); */
+ likelihood(n,mode2, filenum, 0.0 , minplace, maxplace, means,
+ counts, mean, sd, down, outfile,outpsfile,devnull,devnull,
+ a,b,.20,.80,.03);
+ sclose(outfile);
+ sclose(outpsfile);
+ {
+ char clean_command[200];
+
+ strcpy(clean_command, "rm -f ");
+
+/*** The following removes the pir.txt output file. */
+ system(strcat(clean_command,txt_filename));
+
+ } /** DELETE TEMP FILES SO OTHER PEOPLE CAN **/
+ /** WRITE TO SAME FILESNAMES. **/
+ {
+ char clean_command[200];
+
+ strcpy(clean_command, "rm -f ");
+
+ system(strcat(clean_command,PIR_LIKES));
+
+ }
+
+ *m=a[0]; *b1=b[0];
+}
+
+
+
+
+
+/* Look in file .likelihood for likelihood line. If not there, compute */
+/* it from pir. */
+void get_likelihood_line(char likelihoods[MAX_TABLE_NUMBER][MAXLINE], double *m, double *b, double *m_single, double *b_single, char lib[MAXFUNCTNUM], int functnum, char *pir_name, int mode, int number_tables)
+{
+ char *tmp;
+ int libnumber=0;
+ int i, cnt=0;
+ int current_lib;
+ char buf[80];
+
+ FILE *pir;
+ FILE *linefile;
+ FILE *txt_file1, *txt_file2;
+ double score, scs[MAXSEQLEN][POSNUM];
+ Sequence sequence;
+ char seq[MAXSEQLEN], reg[MAXSEQLEN], offsets[MAXSEQLEN];
+ char title[500], code[500];
+ int seqlen;
+ int likefile;
+ int found_like_lib=0; int found_single=0;
+
+ sequence.seq=seq;
+ sequence.reg=reg;
+ sequence.title=title;
+ sequence.code=code;
+
+
+ for (i=0; i<functnum; i++)
+ libnumber |= (int) pow(2,lib[i]);
+ fprintf(stderr,"\nlibnumber= %d",libnumber);
+
+ /* What directory to look for paircoil defauts in. Set environment */
+ /* variable LIKELINE_DIR to be a path to this file. */
+ for (likefile=0; likefile <number_tables; likefile++) {
+ current_lib = -1;
+ found_like_lib =0;
+ found_single =0;
+ fprintf(stderr,"\nLikelihood %d=%s",likefile,likelihoods[likefile]);
+ if (likelihoods[likefile][0] == ',') {
+ if (tmp=getenv("LIKELINE_DIR"))
+ strcpy(likelihoods[likefile],tmp);
+ else {
+ strcpy(likelihoods[likefile],getenv("HOME"));
+ strcat(likelihoods[likefile],"/.likelihood");
+ }
+ }
+
+
+
+ if (linefile=sopen(likelihoods[likefile],"r")) {
+ while ( (fgets(buf,80,linefile)) && !(found_single && found_like_lib)) {
+ if (buf[0] == '0') {
+ sscanf(buf,"%*d %lf %lf", &m_single[likefile], &b_single[likefile]);
+ found_single=1;
+ }
+ else if (!found_like_lib) {
+ sscanf(buf, "%d %lf %lf",¤t_lib, &m[likefile], &b[likefile]);
+ if (current_lib == libnumber) found_like_lib=1;;
+ }
+ }
+ sclose(linefile);
+
+ if (!(found_like_lib && found_single)) {
+ if (!found_single)
+ fprintf(stderr,"\nSingles likelihood function %d unknown.",likefile);
+
+ if (current_lib != libnumber ) /* Wasn't in file so should create it
+ by running on pir. Add this later.*/
+ fprintf(stderr, "\nLikelihood function %d unknown...", likefile);
+
+
+ if (pir_name[0]==',') {
+ fprintf(stderr,"\nUnable to continue. Please input the location of the pir in the .paircoil file.");
+ exit(-1);
+ }
+ else {
+ fprintf(stderr, "\nPlease wait a few minutes while the likelihood is computed.");
+
+ pir = sopen(pir_name, "r");
+
+ if (!found_like_lib) txt_file1= sopen(PAIR_PIR_TXT,"w");
+ if (!found_single) txt_file2 = sopen(SINGLE_PIR_TXT,"w");
+
+ for (cnt=1;getseq2(pir,&sequence); cnt++)
+ if (sequence.seqlen >= MINWINDOW) {
+ seqlen= sequence.seqlen;
+ if (!(cnt %10)) {
+ putc('.',stderr);
+ fflush(stderr);
+ }
+
+ WINDOW=PAIRWINDOW;
+ if (!found_like_lib) {
+ switch_tables(likefile);
+ scseqadj(lib,seq,seqlen,scs,offsets, &score,
+ mode,0);
+ txt_output(sequence,0,txt_file1, scs, NULL,offsets,0);
+ }
+ if (!found_single) {
+ switch_tables(likefile);
+
+ scorestock(seq,seqlen,scs,offsets, &score,
+ mode & PROLINE_FREE_MODE, 0);
+ /* Score stock = singles method using our table. */
+ txt_output(sequence,0,txt_file2, scs, NULL,offsets,0);
+ }
+
+ } /* Ends if sequence > MINWINDOW, and hence the for */
+ /* loop that reads from pir. */
+ sclose(pir);
+
+
+ if (!found_like_lib) { /* Now compute the lib likelihood line. */
+ sclose(txt_file1);
+
+ likelihood_from_txt(PAIR_PIR_TXT,&m[likefile],&b[likefile]);
+
+ linefile=sopen(likelihoods[likefile],"a");
+ fprintf(linefile,"%d %lf %lf\n",libnumber,m[likefile], b[likefile]);
+
+ fprintf(linefile,"\t (lib =");
+ /* Now print human readable lib number. */
+ for (i=0; i<functnum; i++)
+ fprintf(linefile," %d",lib[i]);
+ fprintf(linefile," )\n");
+
+ sclose(linefile);
+ }
+
+ if (!found_single) { /* Now compute the singles method like line */
+ sclose(txt_file2);
+
+ likelihood_from_txt(SINGLE_PIR_TXT,&m_single[likefile],
+ &b_single[likefile]);
+
+ linefile=sopen(likelihoods[likefile],"a");
+ fprintf(linefile,"%d %lf %lf\n",0 ,m_single[likefile],
+ b_single[likefile]);
+
+ fprintf(linefile,"\t (lib =");
+ /* Now print human readable lib number. */
+ fprintf(linefile," Singles Method");
+ fprintf(linefile," )\n");
+
+ sclose(linefile);
+ }
+ } /* Ends else from if (pir_name == ',') */
+
+ } /* Ends section for file not having lib or singles like line */
+
+ } /* Ends section where search likelihoods[likefile] file. */
+ else {
+ fprintf(stderr,"\nCouldn't find likelihood line file %d\n",likefile);
+ exit(-1);
+ }
+ } /* Ends count through tables=linefiles. */
+}
+
+
+
+double stocksc2likelihood (double sc)
+{
+ double gg, gcc;
+
+ sc = exp(sc/28);
+ gg = (sc - NEGMEAN)/NEGSD;
+ gcc = (sc - CCMEAN)/CCSD;
+ gg = exp((-gg * gg) / 2) / NEGSD / SQRT2PI;
+ gcc = exp((-gcc * gcc) / 2) / CCSD / SQRT2PI;
+ return(gcc / (NEGCCRATIO * gg + gcc));
+}
+
+
+
+/* Likelihood line is mx +b */
+double sc2likelihood (double sc, double m, double b)
+{
+ double likelihood;
+
+ likelihood = m*sc + b;
+
+ if (likelihood>1) return(1);
+ else if (likelihood<0) return(0);
+ else return(likelihood);
+}
+
+
+static int current_seq_number=0;
+
+void PairCoilScore (int mode, Sequence sequence, char lib[MAXFUNCTNUM], double m, double b, double *maxscore, int table, char offsets[MAXSEQLEN], double scores[MAXSEQLEN][POSNUM+1], int offset_to_use, int raw_score)
+{
+ int i,j;
+ double sc2scores[MAXSEQLEN][POSNUM];
+
+ int ProlineFreeWindow= mode & PROLINE_FREE_MODE;
+
+ double frac_coiled_this_table;
+ extern double init_class_prob[MAX_CLASS_NUMBER];
+
+ frac_coiled_this_table = init_class_prob[table];
+
+
+ switch_tables(table);
+
+ /* score the sequence using sc2seq */
+ scseqadj(lib,
+ sequence.seq, sequence.seqlen,
+ sc2scores,
+ offsets,
+ maxscore, mode, 0);
+
+
+ if (!raw_score)
+ if (mode & USE_LIKE_LINE)
+ *maxscore= sc2likelihood(*maxscore,m,b);
+ else {
+ *maxscore= frac_coiled_this_table * exp(*maxscore);
+ }
+
+ for (i=0; i<sequence.seqlen; ++i) {
+ /* Do the most probabable register if offset_to_use is -1 or 7. */
+ if ((offset_to_use == -1) || (offset_to_use ==7))
+ if (raw_score) /*** Return rawscore. **/
+ scores[i][7] = sc2scores[i][offsets[i]];
+ else
+ if (mode & USE_LIKE_LINE)
+ scores[i][7] = sc2likelihood(sc2scores[i][offsets[i]],m,b);
+ else {
+ scores[i][7] = frac_coiled_this_table *exp(sc2scores[i][offsets[i]]);
+ if (scores[i][7] >1) scores[i][7]=1;
+ }
+
+ else /* Offset to use is 0-6. */
+ if (raw_score)
+ scores[i][7] = sc2scores[i][offset_to_use];
+ else
+ if (mode & USE_LIKE_LINE)
+ scores[i][7] = sc2likelihood(sc2scores[i][offset_to_use],m,b);
+ else {
+ scores[i][7] = frac_coiled_this_table *
+ exp(sc2scores[i][offset_to_use]);
+ if (scores[i][7] >1) scores[i][7]=1;
+ }
+
+ /* do other registers */
+ for (j=0; j<7; ++j) { /** sc2scores[i][j] is the score for position
+ i+j%POSNUM **/
+ if (raw_score)
+ scores[i][(i+j)%POSNUM] = sc2scores[i][j];
+ else
+ if (mode & USE_LIKE_LINE)
+ scores[i][(i+j)%POSNUM] = sc2likelihood(sc2scores[i][j],m,b);
+ else {
+ scores[i][(i+j)%POSNUM] = frac_coiled_this_table *
+ exp(sc2scores[i][j]);
+ if (scores[i][(i+j)%POSNUM] >1) scores[i][(i+j)%POSNUM]=1;
+ }
+ }
+ }
+
+}
+
+
+
+
+
+/***************************************************************************/
+
+/** Score using the singles method on our tables. **/
+double *SingleCoilScore(int mode, Sequence sequence, double m, double b,
+ double *maxscore, int table,
+ double scores[MAXSEQLEN][POSNUM+1], int offset_to_use,
+ int raw_score)
+{
+ int i,j;
+ double singlescores[MAXSEQLEN][POSNUM];
+ char offsets[MAXSEQLEN];
+
+ int ProlineFreeWindow= mode & PROLINE_FREE_MODE;
+
+ double frac_coiled_this_table;
+
+
+ switch_tables(table);
+
+ /* score the thing */
+ scorestock(
+ sequence.seq, sequence.seqlen,
+ singlescores,
+ offsets,
+ maxscore,ProlineFreeWindow, 0); /** Score using our table. */
+
+
+ if (!raw_score)
+ if (mode & USE_LIKE_LINE)
+ *maxscore= sc2likelihood(*maxscore,m,b);
+ else {
+ /* fprintf(stderr,"\nmaxscore = %lf", *maxscore);
+ */
+ *maxscore= frac_coiled_this_table * exp(*maxscore);
+ /* fprintf(stderr," maxlike = %lf", *maxscore);
+ */
+ }
+
+ for (i=0; i<sequence.seqlen; ++i) {
+ /* Do the most likely register. */
+ if ( (offset_to_use == -1) || (offset_to_use == 7)) /* Use max offset. */
+ if (!raw_score)
+ if (mode & USE_LIKE_LINE)
+ scores[i][7] = sc2likelihood(singlescores[i][offsets[i]],m,b);
+ else {
+ scores[i][7] = frac_coiled_this_table*
+ exp(singlescores[i][offsets[i]]);
+ if (scores[i][7] >1) scores[i][7]=1;
+ }
+ else scores[i][7] = singlescores[i][offsets[i]];
+
+ else /* Use offset_to_use. */
+ if (!raw_score)
+ if (mode & USE_LIKE_LINE)
+ scores[i][7] = sc2likelihood(singlescores[i][offset_to_use],m,b);
+ else {
+ scores[i][7] = frac_coiled_this_table*
+ exp(singlescores[i][offsets[i]]);
+ if (scores[i][7] >1) scores[i][7]=1;
+ }
+ else scores[i][7] = singlescores[i][offset_to_use];
+
+ /* do the other registers */
+ for (j=0; j<7; ++j)
+ if (!raw_score)
+ if (mode & USE_LIKE_LINE)
+ scores[i][(i+j)%POSNUM] = sc2likelihood(singlescores[i][j],m,b);
+ else {
+ scores[i][(i+j)%POSNUM] = frac_coiled_this_table*
+ exp(singlescores[i][j]);
+ if (scores[i][(i+j)%POSNUM] >1) scores[i][(i+j)%POSNUM]=1;
+ }
+ else scores[i][(i+j)%POSNUM] = singlescores[i][j]; /*Raw_score*/
+
+ }
+
+
+}
+/***************************************************************************/
+
+
+
+void MultiCoilDimensionScore (int mode, Sequence sequence, char lib[MAXFUNCTNUM], double *maxscore, int table, char offsets[MAXSEQLEN], double scores[MAXSEQLEN][POSNUM+1], int offset_to_use)
+{
+ int i,j;
+ double sc2scores[MAXSEQLEN][POSNUM];
+
+ int ProlineFreeWindow= mode & PROLINE_FREE_MODE;
+
+
+ switch_tables(table);
+
+ /* score the sequence using sc2seq */
+ scseqadj(lib,
+ sequence.seq, sequence.seqlen,
+ sc2scores,
+ offsets,
+ maxscore, mode, 0);
+
+/*** Do a SingleCoil type score. ****/
+/***
+ scorestock(sequence.seq, sequence.seqlen,
+ sc2scores,offsets,
+ maxscore,ProlineFreeWindow, 0); **/ /** Score using our table. */
+
+
+ for (i=0; i<sequence.seqlen; ++i) {
+ /* Do the most probabable register if offset_to_use is -1 or 7. */
+ if ((offset_to_use == -1) || (offset_to_use ==7))
+ scores[i][7] = sc2scores[i][offsets[i]];
+ else /* Offset to use is 0-6. */
+ scores[i][7] = sc2scores[i][offset_to_use];
+
+ /* do other registers */
+ for (j=0; j<7; ++j) { /** sc2scores[i][j] is the score for position
+ i+j%POSNUM **/
+ scores[i][(i+j)%POSNUM] = sc2scores[i][j];
+ }
+ }
+
+}
+
+
+
+
+
+/***************************************************************************/
+
+/*** Score using the singles method on stock table and likelihood. **/
+void NEWCOILSScore(int mode, Sequence sequence, double *maxscore,
+ double scores[MAXSEQLEN][POSNUM+1],
+ int offset_to_use)
+{
+ int i,j;
+ double stockscores[MAXSEQLEN][POSNUM];
+ char offsets[MAXSEQLEN];
+
+
+ int ProlineFreeWindow= mode & PROLINE_FREE_MODE;
+
+ /* score the thing */
+ scorestock(
+ sequence.seq, sequence.seqlen,
+ stockscores,
+ offsets,
+ maxscore,ProlineFreeWindow, 1); /** Score using STOCK table. */
+
+
+ *maxscore= stocksc2likelihood(*maxscore);
+
+
+ for (i=0; i<sequence.seqlen; ++i) {
+ /* Do the most likely register. */
+ if ( (offset_to_use == -1) || (offset_to_use == 7)) /* Use max offset. */
+ scores[i][7] = stocksc2likelihood(stockscores[i][offsets[i]]);
+ else /* Use offset_to_use. */
+ scores[i][7] = stocksc2likelihood(stockscores[i][offset_to_use]);
+ /* do the other registers */
+ for (j=0; j<7; ++j)
+ scores[i][(i+j)%POSNUM] = stocksc2likelihood(stockscores[i][j]);
+ }
+
+
+}
+
+
+/***********************************************************/
+
+
+
+
+
+
+void ActualCoils (Sequence sequence, double scores[MAXSEQLEN][POSNUM+1],
+ int offset_to_use, int preprocessor)
+{
+ int i,j;
+
+ if (!sequence.reg)
+ fprintf(stderr,"\nInput file contains no coil regions for sequence in ActualCoil scoring method.");
+
+
+ for (i=0; i<sequence.seqlen; i++) {
+ for (j=0; j<POSNUM; j++)
+ if (!sequence.reg) scores[i][j] =0; /* No coils in input file. */
+
+ else if ( (offset_to_use ==7) || /* max offset is pos file reg */
+ ((offset_to_use == -1) && !preprocessor) )
+ if (j== sequence.reg[i])
+ scores[i][j] = 1;
+ else scores[i][j]=0;
+
+ else if ( (offset_to_use == -1) && preprocessor) /* Do all registers */
+ if (sequence.reg[i] != '.')
+ scores[i][j] = 1;
+ else scores[i][j]=0;
+
+ else /* just give likelihood 1 to offset_to_use. */
+ if ( (sequence.reg[i] != '.') && (j == (i+offset_to_use)%POSNUM))
+ scores[i][j] = 1;
+ else scores[i][j]=0;
+
+ if (!sequence.reg) scores[i][7]=0;
+ else if (sequence.reg[i] != '.')
+ scores[i][7]=1;
+ else scores[i][7]=0;
+ }
+
+}
+
+
+
+
+
+
+
+/********************************************************/
+
+/*** This averages score over coil, where a register shift does not start **/
+/*** a new coil if start_at_reg_shift =0. ****/
+/*** This defines whether in a coil by if likelihood is > bound. ***/
+/** Value of 0 for avg_max means do average, 1 means do max. ***/
+/*** For max score, find the maximum scoring residue in the CoilLike ***/
+/*** and take the corresponding StructureScore. ***/
+/*** To compute the max version, average over all residues in the coil **/
+/*** that score the max coil likelihood. ***/
+/*** BE VERY CAREFUL in making changes!!! i.e. need to count down on k **/
+/*** since offset 7 needs to make use of previous offsets, so don't want **/
+/*** them to have been changed to coil score yet. Also, we sometimes ***/
+/*** pass in CoilLike as the StructureLike too, in which case have to be**/
+/*** doubly careful. There is also ONE place where make a reference **/
+/** to offset 7 when dealing with other registers, so had to save an **/
+/** old version of offset 7 in case it changes (when StructureLike **/
+/** is the same as CoilLike. **/
+/*** ACTUALLY, all that should be moot now that we changed it so residue **/
+/*** scores are always saved in tablenum +3, so pass those in and save **/
+/*** coil scores in tablenum. ***/
+void average_score_over_coils3(Sequence seq,
+ double StructureResLike[MAXSEQLEN][POSNUM+1],
+ double StructureCoilLike[MAXSEQLEN][POSNUM+1],
+ double CoilLike[MAXSEQLEN][POSNUM+1],
+ int offset_to_use,
+ char offsets[MAXSEQLEN],
+ int start_at_reg_shift, double bound,
+ int avg_max)
+{
+ int i=0, j,k; /* k is the offset. */
+ int start_coil[POSNUM+1];
+ double total_coil_score[POSNUM+1];
+ double avg_coil_score[POSNUM+1];
+ double max_CoilLike[POSNUM+1];
+ double total_max_StructLike[POSNUM+1];
+ int number_max_residues[POSNUM+1];
+ double max_StructLike[POSNUM+1];
+ int regist;
+
+ double total_coil_length[POSNUM +1]; /* A weighted length. */
+ double CoilLike_MaxOffset[MAXSEQLEN];
+
+
+
+ for (j=0; j<POSNUM+1; j++) {
+ start_coil[j] = -1;
+ total_coil_score[j]=0;
+ avg_coil_score[j]=0;
+ total_coil_length[j]=0;
+ max_CoilLike[j]=0;
+ total_max_StructLike[j]=0;
+ number_max_residues[j]=0;
+ }
+
+ while (i<=seq.seqlen) {
+ if (i<seq.seqlen) CoilLike_MaxOffset[i]=CoilLike[i][7];
+
+ for (k=POSNUM; k>=0; k--) { /** Count down since k=7 needs to use **/
+ /** scores from other offset (i.e see **/
+ /** use of regist_to_use). **/
+ if (k!=7) regist = (i+k)%POSNUM;
+ else regist =7;
+
+ /* Consider <= bound to be out of coil */
+ if ( ( CoilLike[i][regist] <= bound) || (i==seq.seqlen) ||
+ ( (i>0) && (offsets[i] != offsets[i-1]) && start_at_reg_shift &&
+ (CoilLike_MaxOffset[i-1]> bound) ) ) {
+ /** If not in a coil (value of -HUGE_VAL signals not in coil). */
+ if (start_coil[k] != -1) { /* Ended coil at i-1. */
+ avg_coil_score[k] = total_coil_score[k]/total_coil_length[k];
+ max_StructLike[k] = total_max_StructLike[k]/number_max_residues[k];
+ for (j=start_coil[k]; j<i; j++) {
+ if (k!=7)
+ if (avg_max == 0)
+ StructureCoilLike[j][(j+k)%POSNUM]= avg_coil_score[k];
+ else
+ StructureCoilLike[j][(j+k)%POSNUM] = max_StructLike[k];
+ else
+ if (avg_max == 0)
+ StructureCoilLike[j][7]= avg_coil_score[k];
+ else
+ StructureCoilLike[j][7]= max_StructLike[k];
+ }
+ }
+ total_coil_score[k] =0; /***Resets for "not in coil" ***/
+ max_CoilLike[k]=0;
+
+ if (i != seq.seqlen) /* Consider <= bound to be out of coil */
+ if ( CoilLike[i][regist] <= bound) { /* Not in coil anymore. */
+ StructureCoilLike[i][regist]=0; /** Make coil score under bound */
+ /** be 0**/
+
+/*** StructureCoilLike[i][regist]= StructureResLike[i][regist]; **/
+ /** If not in coil above bound make **/
+ /** the coil like the same as the **/
+ /** regist like. **/
+ start_coil[k] = -1;
+ total_coil_length[k] =0;
+ max_CoilLike[k]=0;
+ }
+ else { /* Changing registers. */
+ start_coil[k]= i;
+ if (CoilLike[i][regist] > max_CoilLike[k]) {
+ max_CoilLike[k] = CoilLike[i][regist];
+ total_max_StructLike[k] = StructureResLike[i][regist];
+ number_max_residues[k] =1;
+ }
+ else if (CoilLike[i][regist] == max_CoilLike[k]) {
+ total_max_StructLike[k] += StructureResLike[i][regist];
+ number_max_residues[k]++;
+ }
+
+ if (1) { /** Do weighted_avg **/
+ int reg_to_use=regist;
+ if (regist ==7)
+ reg_to_use= (i+offsets[i])%POSNUM;
+
+ total_coil_length[k]= CoilLike[i][reg_to_use];
+ total_coil_score[k] = StructureResLike[i][regist]*
+ CoilLike[i][reg_to_use];
+/*********
+ if (total_coil_length[k] ==0)
+ fprintf(stderr,"\n0 coil length at %d, reg_to_use %c",
+ i, 'a'+reg_to_use);
+*************/
+ }
+ }
+ } /** Ends if ending a coil. **/
+
+ else { /* If in a coil. */
+ if (CoilLike[i][regist] > max_CoilLike[k]) {
+ max_CoilLike[k] = CoilLike[i][regist];
+ total_max_StructLike[k] = StructureResLike[i][regist];
+ number_max_residues[k] =1;
+ }
+ else if (CoilLike[i][regist] == max_CoilLike[k]) {
+ total_max_StructLike[k] += StructureResLike[i][regist];
+ number_max_residues[k]++;
+ }
+
+ if (1) { /* Weighted average scores. */
+ int reg_to_use=regist;
+ if (regist ==7)
+ reg_to_use= (i+offsets[i])%POSNUM;
+
+ total_coil_score[k] += StructureResLike[i][regist]*
+ CoilLike[i][reg_to_use];
+ total_coil_length[k] += CoilLike[i][reg_to_use];
+ }
+ if (start_coil[k] == -1)
+ start_coil[k] =i;
+ }
+ }
+
+ i++;
+ }
+
+ /* Need the following hack to get the correct register in combo off method.*/
+/************
+ if ( (offset_to_use == 7) && (!start_at_reg_shift))
+ for (i=0; i<seq.seqlen; i++) {
+ for (j=0; j<POSNUM; j++)
+ StructureCoilLike[i][j]= -HUGE_VAL;
+ StructureCoilLike[i][(i+offsets[i]) %POSNUM] = StructureCoilLike[i][7];
+ }
+ else
+*************/
+
+/** Need the following hack to get scores for max offset method. **/
+ if (offset_to_use == -1)
+ for (i=0; i<seq.seqlen; i++)
+ StructureCoilLike[i][7]= StructureCoilLike[i][(i+offsets[i]) %POSNUM];
+
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
git://source.jalview.org / jpred.git / blobdiff