X-Git-Url: http://source.jalview.org/gitweb/?a=blobdiff_plain;f=sources%2Fmulticoil%2Foutput_cuts.c;fp=sources%2Fmulticoil%2Foutput_cuts.c;h=c58ff2a5e6391f4ac35541c94b27552fed96178d;hb=81362e35a140cd040e948b921053e74267f8a6e3;hp=0000000000000000000000000000000000000000;hpb=2cf032f4b987ba747c04159965aed78e3820d942;p=jpred.git diff --git a/sources/multicoil/output_cuts.c b/sources/multicoil/output_cuts.c new file mode 100644 index 0000000..c58ff2a --- /dev/null +++ b/sources/multicoil/output_cuts.c @@ -0,0 +1,1201 @@ +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