changed void finish_log(int mode, double bound, FILE *fin, FILE *flog, FILE *flog_coil_conflicts, FILE *fout_coils, FILE *fout, /* FILE *fout_raw[MAX_TABLE_NUMBER+1], */ double *m, double *b, double *m_single, double *b_single, char lib[MAXFUNCTNUM], char lib_differ[MAXFUNCTNUM], char multi_lib[MAX_TABLE_NUMBER][MAXFUNCTNUM], int main_method, int main_table, int main_preprocessor_method, int *seqcnt, int avg_max, int which_priors, int which_priorp, double prior_freq_single[MAX_TABLE_NUMBER], double prior_freq_pair[MAX_TABLE_NUMBER], int good_turing_fixed_disc, int structural_pos[POSNUM+1]) to void finish_log(int mode, double bound, FILE *fin, FILE *flog, FILE *fout_coils, FILE *fout, double *m, double *b, double *m_single, double *b_single, char lib[MAXFUNCTNUM], char multi_lib[MAX_TABLE_NUMBER][MAXFUNCTNUM], int main_method, int main_table, int main_preprocessor_method, int *seqcnt, int avg_max) changed output_seq(lib,differ_lib, multi_lib,m,b,m_single, b_single, mode,bound,flog,flog_coil_conflicts, fout_coils,fout, avg_max, main_method, main_preprocessor_method, main_table); to output_seq(lib, multi_lib,m,b,m_single, b_single, mode,bound,flog, fout_coils,fout, avg_max, main_method, main_preprocessor_method, main_table); changed NewSeq_nonX (mode, sequence, which_priors, which_priorp, prior_freq_single, prior_freq_pair, good_turing_fixed_disc, structural_pos); to NewSeq_nonX (mode, sequence); cut /*** In TableDiffCoil method and PairDiffAvg method ***/ /*** this computes the distribution of coil scores. ***/ /*** If in pos mode it only counts coils that are ***/ /*** in the pos file. ***/ /*** It is called from output_seq. The results ***/ /*** are printed from finish_log to log file. ***/ void count_positives(int *number_coils_at_most_half, int *number_coils_above_half, int *number_coils_above_3quarters, int *number_coils_below_1quarter, int *number_coils_missed_by_preproc, Sequence seq, double scores[MAXSEQLEN][POSNUM+1], double preproc_like[MAXSEQLEN][POSNUM+1], int mode, double bound, int start_at_reg_shift, FILE *flog, double *minpositive) /* minpositive[0] is min differentiator score on a pos coil. */ /* minpositive[1] is min differentiator score above .5 on pos coil. */ /* minpositive[2] is min preproc score above bound on a pos coil. */ { int i,j; char offsets[MAXSEQLEN]; int not_scored_current_coil =1; int in_pos_coil =0; int printed_to_log =0; double preproc_score_current_coil =0; for (i=0; i0) && (offsets[i-1] != -1) && (offsets[i] != offsets[i-1]) && (preproc_like[i][(i+offsets[i-1]) %POSNUM] == preproc_like[i][7]) ) offsets[i] = offsets[i-1]; } if (mode & POS_MODE) { if (seq.reg) /* Only count coils if in posfile */ for (i=0; i preproc_score_current_coil) preproc_score_current_coil = preproc_like[i][7]; /* Score of */ /* coil is max over coil */ in_pos_coil=1; if (preproc_like[i][7]>bound) { if ( not_scored_current_coil || /* A new coil. */ ( (offsets[i] != offsets[i-1]) && start_at_reg_shift)) { if (scores[i][7] < minpositive[0]) minpositive[0] = scores[i][7]; if ((scores[i][7] < minpositive[1]) && (scores[i][7]>.5)) minpositive[1] = scores[i][7]; if (scores[i][7]>.5) { (*number_coils_above_half)++; if (scores[i][7] > .75) (*number_coils_above_3quarters)++; } else { printed_to_log=1; fprintf(flog,"ALERT: Low score at res %d in %s\n",i,seq.code); (*number_coils_at_most_half)++; if (scores[i][7] < .25) (*number_coils_below_1quarter)++; } not_scored_current_coil =0; } } else { /* Coil in pos file, but not found by preproc. */ if ( start_at_reg_shift && not_scored_current_coil && (seq.reg[i]%POSNUM != ((seq.reg[i-1]+1) %POSNUM)) && (seq.reg[i-1] != '.') ) { /* Register shift in pos */ (*number_coils_missed_by_preproc)++; printed_to_log=1; fprintf(flog,"ALERT: Missed Coil at res %d in %s\n",i,seq.code); not_scored_current_coil =1; } } } else { /* Not a posfile coil. */ if ((in_pos_coil) && (preproc_score_current_coil < minpositive[2]) && (preproc_score_current_coil>bound) ) minpositive[2] = preproc_score_current_coil; if ( (not_scored_current_coil) && (in_pos_coil)) { /* Never found current coil. */ (*number_coils_missed_by_preproc)++; printed_to_log=1; fprintf(flog,"ALERT: Missed Coil at res %d in %s\n",i,seq.code); } not_scored_current_coil =1; in_pos_coil =0; } } if ((in_pos_coil) && (preproc_score_current_coil < minpositive[2]) && (preproc_score_current_coil>bound) ) /* If last pos coilended at */ /* end of seq. */ minpositive[2] = preproc_score_current_coil; if (not_scored_current_coil && in_pos_coil) { /* missed last coil. */ (*number_coils_missed_by_preproc)++; printed_to_log=1; fprintf(flog,"ALERT: Missed Coil at res %d in %s\n",i,seq.code); } if (printed_to_log) fprintf(flog,"\n\n\n"); } else /* Not in pos mode, so score all coils found. */ for (i=0; ibound) { if ( not_scored_current_coil || /* A new coil. */ ( (offsets[i] != offsets[i-1]) && start_at_reg_shift)) { if (scores[i][7] < minpositive[0]) minpositive[0] = scores[i][7]; if ((scores[i][7] < minpositive[1]) && (scores[i][7]>.5)) minpositive[1] = scores[i][7]; if (scores[i][7]>.5) { (*number_coils_above_half)++; if (scores[i][7] > .75) (*number_coils_above_3quarters)++; } else { (*number_coils_at_most_half)++; if (scores[i][7] < .25) (*number_coils_below_1quarter)++; } not_scored_current_coil =0; } } else /* Nothing found by preproc. */ not_scored_current_coil =1; } } /*********************************ROUND 2*********************************/ /** Used to print out for deciding between 2 and 3 stranded. */ void log_output3_prn(int mode, double score, double scs[MAXSEQLEN][POSNUM+1], int seqlen, int cnt, Sequence sequence, double bound, double upper_bound, FILE *flog) { char *title, *code, *reg, *seq; int i, st; double max_pos, max_neg; double the_score; char regist, offset; char offsets[MAXSEQLEN]; title=sequence.title; code=sequence.code; reg=sequence.reg; seq=sequence.seq; /* Print out sequence number, code, title, and maxscore if > bound. */ /* Note that score represents the max-likelihood by paircoil. */ if (score > 0) { /* As long as there was a coil found by paircoil. */ fprintf(flog, "%d \t%s \n\t%s \n\t%6.2lf \n", cnt, code, title, score); if (sequence.reg && (mode & DEBUG_MODE)) /* Print out coil predictions like a pos file. */ print_coil2(mode, seqlen, reg, scs, bound, upper_bound, seq, flog); /* Gives register predicted by */ /* paircoil. */ /* Now put scores for all the regions that score above bound. */ /* Note that if st<0 then we are not currently in a coiled region. */ fprintf(flog, " ["); st = -1; max_pos = -HUGE_VAL; /* max==-HUGE_VAL and st==-1 equivalently mean that */ max_neg = HUGE_VAL; /* aren't currently in a coil. */ for (i = 0; i < seqlen; i++) { the_score = -HUGE_VAL; for (offset=0; offset0) && (offsets[i-1] != -1) && (offsets[i] != offsets[i-1]) && (scs[i][(i+offsets[i-1]) %POSNUM] == scs[i][7]) ) offsets[i] = offsets[i-1]; if (the_score != -HUGE_VAL ) { /* In a coil*/ if (st < 0) st = i; /* Start new coil at i after non-coiled region.*/ else /* else if i starts a new coil by changing register, */ /* then print previous coil and start new coil. */ if (offsets[i] != offsets[i-1]) { if (max_pos > -HUGE_VAL) { /* ended a coil at i-1 to print out. */ if (max_pos != max_neg) fprintf(flog,"%6.2lf,%6.2lf@%4d-%4d:%c; ", max_neg,max_pos,st,i-1, numpos( (offsets[st] + st) %POSNUM)); else fprintf(flog,"%6.2lf@%4d-%4d:%c; ", max_pos,st,i-1, numpos( (offsets[st] + st) %POSNUM)); } max_pos = scs[i][7]; /*Reset for new coil starting at i. */ max_neg = scs[i][7]; st = i; /* Start new coil at i. */ } if (max_pos < scs[i][7]) max_pos = scs[i][7]; if (max_neg > scs[i][7]) max_neg = scs[i][7]; /* New max for current coil. */ } else if (max_pos > -HUGE_VAL) { /* If just ended a coil by scoring low. */ if (max_pos != max_neg) fprintf(flog,"%6.2lf@,%6.2lf%4d-%4d:%c; ",max_neg,max_pos,st,i-1, numpos((offsets[st]+st) %POSNUM) ); else fprintf(flog,"%6.2lf%4d-%4d:%c; ",max_pos,st,i-1, numpos((offsets[st]+st) %POSNUM) ); st = -1; max_pos = -HUGE_VAL; max_neg = HUGE_VAL; } } /* ends count on i. */ if (st >= 0) /** If the end of the seq. ended a coiled region. */ if (max_pos > -HUGE_VAL) { if (max_pos != max_neg) fprintf(flog,"%6.2lf,%6.2lf@%4d-%4d:%c; ",max_neg,max_pos,st,i-1, numpos((offsets[st]+st) %POSNUM)); else fprintf(flog,"%6.2lf@%4d-%4d:%c; ",max_pos,st,i-1, numpos((offsets[st]+st) %POSNUM)); } fprintf(flog, "End: %d]\n\n", seqlen-1); } /* Ends "if (score > bound)" */ } /****************************************************************************/ void tuple_output(Sequence sequence, int mode, FILE *fout, double all_scs[MAX_NUM_SCORE_DIM][MAXSEQLEN][POSNUM], double all_like[MAX_NUM_SCORE_DIM][MAXSEQLEN][POSNUM+1], char offsets[MAX_NUM_SCORE_DIM][MAXSEQLEN], int number_tables, int number_multi_lib, double bound) { int start, end = -1; int j; char *seq= sequence.seq; char *reg= sequence.reg; int seqlen= sequence.seqlen; int table; int posp = mode & POS_MODE; int above_bound_mode = mode & ABOVE_BOUND_MODE; int res_above_bound; extern Sequence sequence; /* extern FILE *debug_out; */ int dist; if (posp && (!reg)) { /*** fprintf(stderr,"\nIgnored a sequence without posfile registers in txt_output."); **/ return; } /****************************************************************************/ void tuple_output(Sequence sequence, int mode, FILE *fout, double all_scs[MAX_NUM_SCORE_DIM][MAXSEQLEN][POSNUM], double all_like[MAX_NUM_SCORE_DIM][MAXSEQLEN][POSNUM+1], char offsets[MAX_NUM_SCORE_DIM][MAXSEQLEN], int number_tables, int number_multi_lib, double bound) { int start, end = -1; int j; char *seq= sequence.seq; char *reg= sequence.reg; int seqlen= sequence.seqlen; int table; int posp = mode & POS_MODE; int above_bound_mode = mode & ABOVE_BOUND_MODE; int res_above_bound; extern Sequence sequence; /* extern FILE *debug_out; */ int dist; if (posp && (!reg)) { /*** fprintf(stderr,"\nIgnored a sequence without posfile registers in txt_output."); **/ return; } while (nondot_advance(&start, &end, seq, reg, seqlen, posp)) { if (end-start < MINWINDOW) { continue; } /* If a coil is shorter than */ /* MINWINDOW, it does not get*/ /* output to the histogram. */ /*** Note -HUGE_VAL means was a proline, should be -HUGE_VAL for all tables */ for (j=start;jbound); if (res_above_bound && (!posp || reg[j] != '.')) { /*************** if ( (all_scs[0][j][offsets[0][j]] < 0) && (all_scs[1][j][offsets[1][j]] < 8) ) fprintf(debug_out,"In sequence %s %s residue %d scored %lf, %lf\n", sequence.code, sequence.title, j, all_scs[0][j][offsets[0][j]],all_scs[1][j][offsets[1][j]]); ****************/ fprintf(fout,"\n"); for (dist=0; dist < number_multi_lib; dist++) for (table =0; tablecode_for_sequence1 start1 end1 */ /* >code_for_sequence2 start2 end2 */ /* .... */ res_compare (FILE *ps, double *sc1, double *sc2, Sequence *sequence, FILE *regionfile) { char c; int i; int reset; char codebuff[MAXLINE]; int rstart, rend, seqlen; double max1=-HUGE_VAL, max2=-HUGE_VAL, rmax1=-HUGE_VAL, rmax2=-HUGE_VAL; seqlen=sequence->seqlen; if (regionfile) { while (1) { c=fgetc(regionfile); while (c!='>' && c!=EOF) c=fgetc(regionfile); /* scan for the sequence code name, and a start and end value */ if(c=='>') { fscanf (regionfile, "%s %d %d", codebuff, &rstart, &rend); if ( !strcmp (codebuff, sequence->code)) break; /** We got the code for the sequence we are working on */ } else { rstart=1; rend=seqlen; break; /** No more codes in the graph file to check, so just */ /* break and graph the whole sequence. */ } } --rstart; } else { /* Just graph the whole seq if no region file. */ rstart = 0; rend = seqlen; } for (i=rstart; i < rend; ++i) if (rmax1code, sequence->title); fputs (" /sy 125 def /ry 100 def /ty 150 def \n", ps); fputs (" /sx 560 def /ux 650 def /tx 500 def\n", ps); fputs ( "tx sy moveto(sequence) rshow tx ry moveto(region) rshow\n", ps); fputs ("sx ty moveto(Stock) cshow ux ty moveto(us) cshow\n", ps); fprintf (ps, "sx sy moveto(%lf) cshow ux sy moveto(%lf) cshow\n", max1, max2); fprintf (ps, "sx ry moveto(%lf) cshow ux ry moveto(%lf) cshow\n", rmax1, rmax2); fprintf (ps, "315 sy moveto (region %d-%d) show\n", rstart+1, rend ); --rend; --rstart; } else { /* Just graph the whole seq if no region file. */ rstart = 0; rend = seqlen; } for (i=rstart; i < rend; ++i) if (rmax1code, sequence->title); fputs (" /sy 125 def /ry 100 def /ty 150 def \n", ps); fputs (" /sx 560 def /ux 650 def /tx 500 def\n", ps); fputs ( "tx sy moveto(sequence) rshow tx ry moveto(region) rshow\n", ps); fputs ("sx ty moveto(Stock) cshow ux ty moveto(us) cshow\n", ps); fprintf (ps, "sx sy moveto(%lf) cshow ux sy moveto(%lf) cshow\n", max1, max2); fprintf (ps, "sx ry moveto(%lf) cshow ux ry moveto(%lf) cshow\n", rmax1, rmax2); fprintf (ps, "315 sy moveto (region %d-%d) show\n", rstart+1, rend ); --rend; --rstart; } else { /* Just graph the whole seq if no region file. */ rstart = 0; rend = seqlen; } for (i=rstart; i < rend; ++i) if (rmax1code, sequence->title); fputs (" /sy 125 def /ry 100 def /ty 150 def \n", ps); fputs (" /sx 560 def /ux 650 def /tx 500 def\n", ps); fputs ( "tx sy moveto(sequence) rshow tx ry moveto(region) rshow\n", ps); fputs ("sx ty moveto(Stock) cshow ux ty moveto(us) cshow\n", ps); fprintf (ps, "sx sy moveto(%lf) cshow ux sy moveto(%lf) cshow\n", max1, max2); fprintf (ps, "sx ry moveto(%lf) cshow ux ry moveto(%lf) cshow\n", rmax1, rmax2); fprintf (ps, "315 sy moveto (region %d-%d) show\n", rstart+1, rend ); --rend; fprintf(ps, "/vx %lf def /vy 5 def\n",(11*72-100)/(double) sequence->seqlen); fputs("/adjust {exch vx mul exch vy mul} def\n",ps); fputs("/-Inf -150 def\n",ps); fputs("/m {adjust moveto vx 0 rlineto} def\n",ps); fputs("/l {adjust lineto vx 0 rlineto} def\n",ps); fputs("/min {1 index 1 index lt {pop} {exch pop} ifelse} def\n",ps); fputs("/bar {exch vx mul 0 moveto 0 exch vy mul rlineto stroke} def\n",ps); fputs("/foo {2 index 2 index 2 index\ 1 index 1 index gt\ {0.3 0.7 1 setrgbcolor pop bar} {0.5 1 0 setrgbcolor exch pop bar} ifelse\n min 1 0 0 setrgbcolor bar} def\n",ps); /* Add axes to PS output */ fputs("/Times-Roman findfont 10 scalefont setfont\n",ps); fputs("1 setlinewidth 0.7 setgray 0 15 adjust moveto 0 -50 adjust lineto stroke \n",ps); for(i=-50; i<=15; i+= 5) { fprintf(ps,"0 %d adjust moveto -5 0 rmoveto 10 0 rlineto stroke \n",i); fprintf(ps,"0 %d adjust moveto \n", i); fprintf(ps,"7 0 rmoveto (%d) show \n",i); } fprintf(ps,"1 setlinewidth 0 setgray %d 15 adjust moveto %d -50 adjust lineto stroke \n",sequence->seqlen, sequence->seqlen); for(i=-50; i<=15; i+= 5) { fprintf(ps,"%d %d adjust moveto -5 0 rmoveto 10 0 rlineto stroke \n",sequence->seqlen, i); fprintf(ps, "%d %d adjust moveto \n",sequence->seqlen, i); fprintf(ps,"-7 0 rmoveto (%5.2f) rshow \n",exp(i/(double)WINDOW)); } fputs("vx 2 div setlinewidth 0.5 setgray\n",ps); fputs("2 setlinewidth \n",ps); fputs("0 setgray [2] 0 setdash \n",ps); for(reset=1 ,i = 0; i seqlen; ++i){ int inf; inf=(sc1 [i] == -HUGE_VAL); if(!inf) fprintf (ps, "%d %lf %s\n" , i, sc1 [i], reset?"m":"l"); reset=inf; } fputs("stroke \n",ps); fputs("0.7 setgray [] 0 setdash\n",ps); for(reset=1 ,i = 0; i seqlen; ++i){ int inf; inf=(sc2 [i] == -HUGE_VAL); if(!inf) fprintf (ps, "%d %lf %s\n" , i, sc2 [i], reset?"m":"l"); reset=inf; } fputs("stroke \n",ps); fprintf(ps,"0 setlinewidth 0 setgray 0 0 adjust moveto %d 0 adjust lineto stroke \n",sequence->seqlen); fprintf(ps,"0 setlinewidth 0 setgray 0 10 adjust moveto %d 10 adjust lineto stroke \n",sequence->seqlen); fprintf(ps,"0 setlinewidth 0 setgray 0 -10 adjust moveto %d -10 adjust lineto stroke \n",sequence->seqlen); for (i=49; i seqlen; i+=50) fprintf(ps, "%d 0 adjust moveto vx 2 div -5 rmoveto 0 10 rlineto stroke\n", i); fprintf(ps, "%d -15 adjust moveto vx 2 div 0 rmoveto 0 30 vy mul rlineto stroke\n", rstart); fprintf(ps, "%d -15 adjust moveto vx 2 div 0 rmoveto 0 30 vy mul rlineto stroke\n", rend); fputs("showpage\n\n", ps); } /**************************************************************************/ /************The following function is intended to give a user friendly ***/ /************output of where the coil predictors and coil distinguishers **/ /************come up with different predictions using different tables. **/ /************If all_scores=preproc_scores, then we are just concerned **/ /************with the coil predictions. Otherwise we want the 2-3 **/ /************stranded predictions too. **/ /*** Note that for 2-3 stranded distinguishing, this is only called if **/ /*** doing coil scores, not residue scores. **/ /*** Note right now in the printing we only deal with the max scoring **/ /*** offset. Maybe make an interactive thing to deal with others? **/ /*** NOTE that we use MAX_TABLE_NUMBER +1 for most things, since we use **/ /*** the number_tables +1 entry for pos file info if in pos mode. **/ /** If avg_or_max is 0 then average residue scores over coil. ***/ /** Otherwise do max residue score for the coil score. ***/ void coil_conflicts(int mode, double all_scores[MAX_TABLE_NUMBER][MAXSEQLEN][POSNUM+1], double all_preproc_like[MAX_TABLE_NUMBER][MAXSEQLEN][POSNUM+1], Sequence sequence, int seqnum, double bound, double bound_preprocessor, int number_tables, FILE *flog_coil_conflicts, int avg_or_max) { double all_preproc_coil_scores[MAX_TABLE_NUMBER][MAXSEQLEN][POSNUM+1]; /* This will give each coil found the avg_max_score over that coil. */ int i,j, tablenum, offset, row; char offsets[MAX_TABLE_NUMBER +1 ][MAXSEQLEN]; int coil_starts[MAX_TABLE_NUMBER +1][MAXNUMCOILS][POSNUM +1]; int coil_ends[MAX_TABLE_NUMBER +1][MAXNUMCOILS][POSNUM +1]; /* coil_...[table][i][offset] holds the start and end of the ith coil */ /* in offset. The end of the coil list is marked by a -1. */ int current_coil_number[MAX_TABLE_NUMBER+1][POSNUM+1]; int number_coils[MAX_TABLE_NUMBER+1][POSNUM+1]; int current_coil; int label_for_coil[MAX_TABLE_NUMBER+1][MAXNUMCOILS]; /* Gives each coil in each table a number, so overlapping */ /* coils have same range of numbers. */ int max_coil_number; int coil_above_bound =0; int label_in_row=0; int printed_blank_line; int number_tables_here; if ( (mode & POS_MODE) && sequence.reg) number_tables_here= number_tables+1; /* Locally add 1 to number_tables */ /* since number_tables entry will */ /* contain pos file info. */ else number_tables_here = number_tables; for (tablenum=0; tablenum0) coil_above_bound =1; if (coil_above_bound) fprintf(flog_coil_conflicts, "\n\n%d \t%s \n\t%s", seqnum, sequence.code, sequence.title); /* if (all_scores == all_preproc_like)*/ /* Just concerned with finding coils */ /*** Put 3 coils per row. ***/ for (row =0; 3*row <= max_coil_number; row++) { printed_blank_line = 0; for (tablenum=0; tablenum0) && (!printed_blank_line)) { fprintf(flog_coil_conflicts,"\n"); /* To print pretty. */ printed_blank_line=1; } if (tablenum != number_tables) /* Regular table. */ fprintf(flog_coil_conflicts,"\nTab%d", tablenum); else /* Posfile coils. */ fprintf(flog_coil_conflicts,"\nPos "); /** Print score@start-end:offset,reg **/ for (i=0; (i<3) && (current_coil < number_coils[tablenum][7]) ; i++) /* 3 coils per row...coil label is 3*row+i */ if ( label_for_coil[tablenum][current_coil] == 3*row + i) { if (tablenum != number_tables) { /* Just a regular table. */ if (all_scores == all_preproc_like) fprintf(flog_coil_conflicts," %5.2lf", all_preproc_coil_scores[tablenum] [coil_starts[tablenum][current_coil][7]] [7]); else /* Both preproc. and coil_diff score */ fprintf(flog_coil_conflicts," %5.2lf,%4.2lf", all_preproc_coil_scores[tablenum] [coil_starts[tablenum][current_coil][7]] [7], all_scores[tablenum] [coil_starts[tablenum][current_coil][7]] [7]); /* Now print coil location, etc. */ fprintf(flog_coil_conflicts,"@%4d-%4d:%d,%c", coil_starts[tablenum][current_coil][7], coil_ends[tablenum][current_coil][7], offsets[tablenum][coil_starts[tablenum][current_coil][7]], 'a' + (( coil_starts[tablenum][current_coil][7] + offsets[tablenum] [coil_starts[tablenum][current_coil][7]]) %POSNUM) ); } else /* Posfile coil so don't print score. */ fprintf(flog_coil_conflicts," %4d-%4d:%d,%c", coil_starts[tablenum][current_coil][7], coil_ends[tablenum][current_coil][7], offsets[tablenum][coil_starts[tablenum][current_coil][7]], 'a' + (( coil_starts[tablenum][current_coil][7] + offsets[tablenum] [coil_starts[tablenum][current_coil][7]]) %POSNUM) ); current_coil_number[tablenum][7]++; current_coil++; } else /** The current_coil does not have current_label. **/ /** Made each coil take (19) 24 spaces + 1 at start. **/ fprintf(flog_coil_conflicts, " "); } label_in_row =0; } } } /**************************************************************************/ /*** This will compute the max scoring coils and their offsets... AND **/ /*** compute the scores for all offsets..... *****/ /*** Note that coils start and end at reg shifts only for max offset. */ /*** offset -1 indicates non-coiled region. */ /*** Also computes the coil_starts[table][i][offset] for the i-th coil **/ /*** found in offset by the table. Also computes the coil ends. **/ /** If in pos mode than store info in slot for number_tables **/ /** If avg_or_max is 0 then average residue scores over coil. ***/ /** Otherwise do max residue score for the coil score. ***/ get_coil_score(Sequence sequence, double bound_preprocessor, double all_preproc_like[MAX_TABLE_NUMBER][MAXSEQLEN][POSNUM+1], double all_preproc_coil_scores[MAX_TABLE_NUMBER][MAXSEQLEN][POSNUM+1], char offsets[MAX_TABLE_NUMBER +1][MAXSEQLEN], int number_tables, int coil_starts[MAX_TABLE_NUMBER+1][MAXNUMCOILS][POSNUM +1], int coil_ends[MAX_TABLE_NUMBER+1][MAXNUMCOILS][POSNUM +1], int number_coils[MAX_TABLE_NUMBER+1][POSNUM+1], int mode, int avg_or_max, double all_scores[MAX_TABLE_NUMBER][MAXSEQLEN][POSNUM+1], double bound) { int i, j, tablenum, offset; double max_coil_score[POSNUM+1]; double total_coil_score[POSNUM+1]; int start_coil[POSNUM+1]; int current_coil[POSNUM +1]; int first_error=1; int started_coil = 0; int pos_coil_number = 0; extern int offset_to_use; int offset_to_use_here; if (offset_to_use ==-1) offset_to_use_here = 7; else offset_to_use_here = offset_to_use; /*********************************************************************/ if ((mode & POS_MODE) && sequence.reg) { /** Compute pos coil info in **/ /** index number_tables +1. **/ for (i=0; i<=sequence.seqlen; i++) { if ( (i0) && (offsets[tablenum][i-1] != -1) && (offsets[tablenum][i] != offsets[tablenum][i-1]) && (all_preproc_like[tablenum][i] [(i+offsets[tablenum][i-1]) %POSNUM] == all_preproc_like[tablenum][i][7]) ) offsets[tablenum][i] = offsets[tablenum][i-1]; } /******************/ /*** Find coils for given offset (7 is max offset). ***/ for (offset=0; offset<=POSNUM; offset++) { if ( (all_preproc_like[tablenum][i][offset] > bound_preprocessor) && ( (all_scores == all_preproc_like) || (all_scores[tablenum][i][offset] > bound) ) ) { /* In coil that meets both bounds, or only has to meet preproc bound. */ /*** For offset 7, max register change so score prev coil. ***/ if ( (start_coil[7] != -1) && (offset == 7) && (i>0) && (offsets[tablenum][i] != offsets[tablenum][i-1])) { for (j = start_coil[7]; j MAXNUMCOILS) { if (first_error) { fprintf(stderr, "\n\nERROR: Found more than max number of coils. Found %d coils in seq %s.",current_coil[7],sequence.title); first_error =0; } } else { coil_starts[tablenum][current_coil[7]][7] = start_coil[7]; coil_ends[tablenum][current_coil[7]][7] = i-1; number_coils[tablenum][7]++; current_coil[7]++; } max_coil_score[7]=0; total_coil_score[7]=0; start_coil[7]=-1; } if (start_coil[offset] == -1) /* Starting a coil. */ start_coil[offset] = i; if (all_preproc_like[tablenum][i][offset] >max_coil_score[offset]) max_coil_score[offset] = all_preproc_like[tablenum][i][offset]; total_coil_score[offset] += all_preproc_like[tablenum][i][offset]; } else { /** Not in a coil. */ all_preproc_coil_scores[tablenum][i][offset]=0; if (start_coil[offset] != -1) { /* ended a coil, so score it. */ for (j = start_coil[offset]; j MAXNUMCOILS) { if (first_error && (offset == offset_to_use_here)) { fprintf(stderr, "\n\nERROR: Found more than max number of coils in offset %d. Found %d coils in %s using table %d",current_coil[offset],sequence.title,tablenum, offset_to_use); first_error=0; } } else { coil_starts[tablenum][current_coil[offset]][offset] = start_coil[offset]; coil_ends[tablenum][current_coil[offset]][offset] = i-1; number_coils[tablenum][offset]++; current_coil[offset]++; } max_coil_score[offset]=0; total_coil_score[offset]=0; start_coil[offset]=-1; } } } /* Count on offset. */ } /* Count on i (position in sequence) */ /************** Now deal with coils ended by end of sequence. *****/ for (offset =0; offset<=POSNUM; offset++) if (start_coil[offset] != -1) { /* ended a coil, with end of seq. */ for (j = start_coil[offset]; j MAXNUMCOILS) { if (first_error) { fprintf(stderr,"\n\nERROR: Found more than max number of coils."); first_error =0; } } else { coil_starts[tablenum][current_coil[offset]][offset] = start_coil[offset]; coil_ends[tablenum][current_coil[offset]][offset] = i-1; number_coils[tablenum][offset]++; current_coil[offset]++; } } /* Ends if started a coil in offset. Also ends count on offset. */ /***********************************************************************/ for (offset= 0; offset<=POSNUM; offset++) { /* Indicate end of coil list with -1 */ coil_starts[tablenum][current_coil[offset]][offset] = -1; coil_ends[tablenum][current_coil[offset]][offset] = -1; current_coil[offset]++; } } /* Count on table. */ } int get_coil_order(int coil_starts[MAX_TABLE_NUMBER+1][MAXNUMCOILS][POSNUM +1], int coil_ends[MAX_TABLE_NUMBER+1][MAXNUMCOILS][POSNUM +1], int number_for_coil[MAX_TABLE_NUMBER+1][MAXNUMCOILS], int number_tables, int number_coils[MAX_TABLE_NUMBER+1][POSNUM+1]) { int current_coil_number=0; int current_coil[MAX_TABLE_NUMBER]; int tablenum; int min_start_current=HUGE_VAL; int max_end_current= -HUGE_VAL; int last_end = -HUGE_VAL; int coil_table; int found_coil =1; for (tablenum=0; tablenum < number_tables; tablenum++) current_coil[tablenum]=0; /*** Find the minimum starting unlabeled coil. For coils with same start **/ /** Take the max length one. For a coil completely covered by the prev.**/ /** coil, give it the same number. */ while (found_coil) { found_coil =0; /* Need to verify are still coils to continue. */ min_start_current=HUGE_VAL; max_end_current= -HUGE_VAL; for (tablenum=0; tablenum< number_tables; tablenum++) { if (current_coil[tablenum] < number_coils[tablenum][7]) { found_coil=1; if (coil_starts[tablenum][current_coil[tablenum]][7] == min_start_current) if (coil_ends[tablenum][current_coil[tablenum]][7] > max_end_current) { coil_table= tablenum; max_end_current= coil_ends[tablenum][current_coil[tablenum]][7]; } if (coil_starts[tablenum][current_coil[tablenum]][7] < min_start_current) { min_start_current=coil_starts[tablenum][current_coil[tablenum]][7]; coil_table= tablenum; max_end_current= coil_ends[tablenum][current_coil[tablenum]][7]; } } } /* Ends count on tablenum */ if (found_coil) { if (max_end_current <= last_end) { number_for_coil[coil_table][current_coil[coil_table]] = current_coil_number -1; /* fprintf(stderr,"\nFound coil in table %d. Coil number is %d", coil_table, current_coil_number-1); fprintf(stderr,"\n start= %d, end=%d", min_start_current, max_end_current); */ } else { number_for_coil[coil_table][current_coil[coil_table]]= current_coil_number; last_end = max_end_current; /* fprintf(stderr,"\nFound coil2 in table %d. Coil number is %d", coil_table, current_coil_number); fprintf(stderr,"\n start= %d, end=%d", min_start_current, max_end_current); */ current_coil_number++; } current_coil[coil_table]++; } /* If found coil*/ } /* Ends while loop. */ return(current_coil_number -1); } Changed: void tuple_output2(Sequence sequence, int mode, FILE *fout, double all_scs[MAX_NUM_SCORE_DIM][MAXSEQLEN][POSNUM+1], int number_tables, int number_multi_lib[MAX_TABLE_NUMBER], double bound, int number_score_dim) to void tuple_output2(Sequence sequence, int mode, FILE *fout, double all_scs[MAX_NUM_SCORE_DIM][MAXSEQLEN][POSNUM+1], int number_tables, double bound,int num_score_dim) and cut int num_score_dim; if (number_score_dim) num_score_dim = number_score_dim; else { num_score_dim =0; for (table=0; table