1 static char rcsid[] = "$Id: iit_get.c 222393 2020-04-13 05:43:57Z twu $";
2 #ifdef HAVE_CONFIG_H
3 #include <config.h>
4 #endif
5
6 #include <stdio.h>
7 #include <stdlib.h>
8 #ifdef HAVE_UNISTD_H
9 #include <unistd.h> /* For getopt */
10 #endif
11 #include <string.h>
12 #include <strings.h> /* For index, rindex */
13 #include <ctype.h>
14 #include <math.h> /* For log */
15 #include "bool.h"
16 #include "mem.h"
17 #include "access.h"
18 #include "intlist.h"
19 #include "univinterval.h"
20 #include "interval.h"
21 #include "iit-read-univ.h"
22 #include "iit-read.h"
23 #include "complement.h"
24 #include "fopen.h"
25 #include "getline.h"
26 #include "parserange.h"
27 #include "getopt.h"
28
29
30 #ifdef DEBUG
31 #define debug(x) x
32 #else
33 #define debug(x)
34 #endif
35
36 /************************************************************************
37 * Program options
38 ************************************************************************/
39
40 static bool value_matches_p = false;
41 static bool user_lowvalue_p = false;
42 static bool user_highvalue_p = false;
43 static double lowval;
44 static double highval;
45
46 static char *fieldstring = NULL;
47 static int fieldint = -1;
48 static bool annotationonlyp = false;
49 static bool exactp = false;
50 static bool sortp = false;
51 static bool signedp = true;
52 static int nflanking = 0;
53 static bool centerp = false;
54 static int centerlength = 0;
55 static bool centeruc = false;
56 static bool runlengthp = false;
57 static bool tallyp = false;
58 static bool zeroesp = false;
59 static bool statsp = false;
60 static bool force_label_p = false;
61 static bool force_coords_p = false;
62
63
64 static struct option long_options[] = {
65 /* Input options */
66 {"lowval", required_argument, 0, 'a'}, /* value_matches_p, user_lowvalue_p, lowval */
67 {"highval", required_argument, 0, 'b'}, /* value_matches_p, user_highvalue_p, highval */
68 {"field", required_argument, 0, 'f'}, /* fieldstring */
69 {"annotonly", no_argument, 0, 'A'}, /* annotationonlyp */
70 {"exact", no_argument, 0, 0}, /* exactp */
71 {"sort", no_argument, 0, 'S'}, /* sortp */
72 {"unsigned", no_argument, 0, 'U'}, /* signedp */
73 {"flanking", required_argument, 0, 'u'}, /* nflanking */
74 {"center", required_argument, 0, 'c'}, /* centerp, centerlength */
75 {"centeruc", no_argument, 0, 'H'}, /* centeruc */
76 {"runlength", no_argument, 0, 'R'}, /* runlengthp */
77 {"tally", no_argument, 0, 'T'}, /* tallyp */
78 {"zeroes", no_argument, 0, 'Z'}, /* zeroesp */
79 {"stats", no_argument, 0, 'N'}, /* statsp */
80 {"label", no_argument, 0, 'L'}, /* force_label_p */
81 {"coords", no_argument, 0, 'C'}, /* force_coords_p */
82
83 /* Help options */
84 {"version", no_argument, 0, 'V'}, /* print_program_version */
85 {"help", no_argument, 0, '?'}, /* print_program_usage */
86 {0, 0, 0, 0}
87 };
88
89 static void
print_program_version()90 print_program_version () {
91 fprintf(stdout,"\n");
92 fprintf(stdout,"iit_get: retrieval utility for Interval Index Trees\n");
93 fprintf(stdout,"Part of GMAP package, version %s\n",PACKAGE_VERSION);
94 fprintf(stdout,"Thomas D. Wu, Genentech, Inc.\n");
95 fprintf(stdout,"Contact: twu@gene.com\n");
96 fprintf(stdout,"\n");
97 return;
98 }
99
100 static void
print_program_usage()101 print_program_usage () {
102 fprintf(stdout,"\
103 Usage: iit_get [OPTIONS...] iitfile query\n\
104 iit_get [OPTIONS...] iitfile query types...\n\
105 iit_get [OPTIONS...] iitfile\n\
106 \n\
107 where query is one of the following forms:\n\
108 \n\
109 chr:start..end\n\
110 chr:start\n\
111 chr:\n\
112 start..end\n\
113 start\n\
114 label\n\
115 \n\
116 Options\n\
117 -f, --field=STRING Show given field part of the annotation\n\
118 -L, --label Interpret query or queries as a label, even if it is numeric\n\
119 -C, --coords Interpret query or queries as coords\n\
120 -A, --annotonly Show annotation lines only (no headers)\n\
121 -S, --sort Sort results by coordinates\n\
122 -U, --unsigned Print all intervals as low..high, even those entered as reverse (high < low)\n\
123 -u, --flanking=INT Show flanking segments on left and right\n\
124 \n\
125 Options for specific IIT formats\n\
126 -a, --lowval=DOUBLE Low bound on a values IIT (default -Inf)\n\
127 -b, --highval=DOUBLE High bound on a values IIT (default +Inf)\n\
128 -c, --center=INT Align reads so given position is centered at given column\n\
129 -H, --centeruc Report only reads with upper-case letter at given position\n\
130 -R, --runlength Report runlength IIT file in tally format\n\
131 -T, --tally Report tally IIT file in tally format\n\
132 -Z, --zeroes Include zeroes in tally format\n\
133 -N, --stats Statistics (count, npositions, mean) of tally format\n\
134 \n\
135 -V, --version Show version\n\
136 -?, --help Show this help message\n\
137 \n\
138 \n\
139 The iit_get program retrieves segments from an iit file that overlap a\n\
140 given coordinate or pair of coordinates. Retrieval is done in\n\
141 logarithmic time (for more details on IIT files, see Wu and Watanabe,\n\
142 Bioinformatics 21:1859-1875, 2005). The start coordinate should be\n\
143 less than or equal to the end coordinate. If they are not, the program\n\
144 will reverse them for you. If only a single coordinate is provided,\n\
145 this is equivalent to providing the same number for the start and end\n\
146 coordinate.\n\
147 \n\
148 The given iit file may contain types (which can be displayed by using\n\
149 the -T flag of the iit_dump program). These types may be used to\n\
150 filter the output of the iit file. Multiple types may be specified,\n\
151 which indicates a disjunctive query, such that iit_get returns entries\n\
152 that match any one of the given tags.\n\
153 \n\
154 If no query is provided on the command line, the program will expect\n\
155 one or more queries from stdin, one per line.\n\
156 \n\
157 See also: iit_store, iit_dump\n\
158 ");
159 return;
160 }
161
162
163 #ifdef __STRICT_ANSI__
164 int getopt (int argc, char *const argv[], const char *optstring);
165 #endif
166
167 static char complCode[128] = COMPLEMENT_LC;
168
169 static void
print_forward(char * sequence,int start,int end)170 print_forward (char *sequence, int start, int end) {
171 int i;
172
173 for (i = start; i <= end; i++) {
174 printf("%c",sequence[i]);
175 }
176 return;
177 }
178
179 static void
print_complement(char * sequence,int end,int start)180 print_complement (char *sequence, int end, int start) {
181 int i;
182
183 for (i = end; i >= start; i--) {
184 printf("%c",complCode[(int) sequence[i]]);
185 }
186 return;
187 }
188
189 static void
print_spaces(int n)190 print_spaces (int n) {
191 while (--n >= 0) {
192 printf(" ");
193 }
194 return;
195 }
196
197
198 /* Need to store just the part of the query specified (e.g., 1..10) */
199 static void
print_interval_centered(char * divstring,Chrpos_T coordstart,int index,IIT_T iit,int fieldint)200 print_interval_centered (char *divstring, Chrpos_T coordstart, int index, IIT_T iit, int fieldint) {
201 Interval_T interval;
202 char *label, *annotation, *restofheader, centerchar;
203 bool allocp;
204 int annotlength, left, centerpos;
205
206 if (fieldint < 0) {
207 annotation = IIT_annotation(&restofheader,iit,index,&allocp);
208 if (allocp == true) {
209 FREE(restofheader);
210 }
211 } else {
212 annotation = IIT_fieldvalue(iit,index,fieldint);
213 allocp = true;
214 }
215 annotlength = strlen(annotation);
216 if (annotation[annotlength-1] == '\n') {
217 annotlength--;
218 }
219
220 interval = IIT_interval(iit,index);
221 left = coordstart - Interval_low(interval); /* + length(query) - queryend */
222 if (Interval_sign(interval) < 0) {
223 centerpos = annotlength-left-1;
224 } else {
225 centerpos = left;
226 }
227 centerchar = annotation[centerpos];
228
229 if (centeruc == true && islower(centerchar)) {
230 if (fieldint >= 0 && allocp == true) {
231 FREE(annotation);
232 }
233 } else {
234 print_spaces(centerlength-left);
235 if (Interval_sign(interval) < 0) {
236 print_complement(annotation,annotlength-1,centerpos+1);
237 printf("[%c]",complCode[(int) centerchar]);
238 print_complement(annotation,centerpos-1,0);
239 } else {
240 print_forward(annotation,0,centerpos-1);
241 printf("[%c]",centerchar);
242 print_forward(annotation,centerpos+1,annotlength-1);
243 }
244 print_spaces(centerlength+left-annotlength);
245 if (fieldint >= 0 && allocp == true) {
246 FREE(annotation);
247 }
248
249 printf("\t");
250 if (Interval_type(interval) > 0) {
251 printf("%s\t",IIT_typestring(iit,Interval_type(interval)));
252 }
253
254 if (divstring != NULL) {
255 if (Interval_sign(interval) < 0) {
256 printf("-%s:",divstring);
257 } else {
258 printf("+%s:",divstring);
259 }
260 }
261
262 if (signedp == false) {
263 printf("%u..%u",Interval_low(interval),Interval_high(interval));
264 } else if (Interval_sign(interval) < 0) {
265 printf("%u..%u",Interval_high(interval),Interval_low(interval));
266 } else {
267 printf("%u..%u",Interval_low(interval),Interval_high(interval));
268 }
269 printf("\t");
270
271 label = IIT_label(iit,index,&allocp);
272 printf("%s",label);
273 if (allocp == true) {
274 FREE(label);
275 }
276 printf("\n");
277 }
278
279 return;
280 }
281
282 static char *
get_total_tally(long int * tally,char * ptr)283 get_total_tally (long int *tally, char *ptr) {
284 int n;
285 char *end;
286
287 if ((end = index(ptr,'\n')) == NULL) {
288 fprintf(stderr,"Premature end of line %s\n",ptr);
289 return 0;
290 }
291 /* fprintf(stderr,"Getting tally for %.*s\n",end-ptr,ptr); */
292
293 while (ptr < end) {
294 while (ptr < end && !isdigit((int) *ptr)) {
295 ptr++;
296 }
297 if (ptr < end) {
298 sscanf(ptr,"%d",&n);
299 (*tally) += n;
300 while (ptr < end && !isspace(*ptr)) {
301 ptr++;
302 }
303 while (ptr < end && isspace(*ptr)) {
304 ptr++;
305 }
306 }
307 }
308
309 return ptr;
310 }
311
312
313 static void
print_line(char * ptr)314 print_line (char *ptr) {
315 while (*ptr != '\0' && *ptr != '\n') {
316 printf("%c",*ptr);
317 ptr++;
318 }
319 return;
320 }
321
322
323
324 /* Need to store just the part of the query specified (e.g., 1..10) */
325 static long int
print_interval_tally(Chrpos_T * lastcoord,char * divstring,Chrpos_T coordstart,Chrpos_T coordend,int indexi,IIT_T iit,bool zeroesp)326 print_interval_tally (Chrpos_T *lastcoord, char *divstring, Chrpos_T coordstart, Chrpos_T coordend,
327 int indexi, IIT_T iit, bool zeroesp) {
328 long int total = 0, subtotal;
329 Interval_T interval;
330 char *annotation, *restofheader, *ptr, *nextptr;
331 bool allocp;
332 Chrpos_T chrpos, intervalend;
333
334 annotation = IIT_annotation(&restofheader,iit,indexi,&allocp);
335
336 interval = IIT_interval(iit,indexi);
337 chrpos = Interval_low(interval);
338 intervalend = Interval_high(interval);
339
340 ptr = annotation;
341
342 if (zeroesp == true) {
343 while (*lastcoord < chrpos) {
344 if (statsp == false) {
345 printf("%s\t%u\t%d\t\n",divstring,*lastcoord,0);
346 }
347 (*lastcoord)++;
348 }
349 }
350
351 while (chrpos < coordstart) {
352 if ((ptr = index(ptr,'\n')) == NULL) {
353 fprintf(stderr,"Premature end of tally from %u to %u\n",
354 Interval_low(interval),Interval_high(interval));
355 return total;
356 } else {
357 ptr++;
358 }
359 chrpos++;
360 }
361
362 while (chrpos <= intervalend && chrpos <= coordend) {
363 subtotal = 0;
364 nextptr = get_total_tally(&subtotal,ptr);
365 if (subtotal > 0 || zeroesp == true) {
366 if (statsp == false) {
367 printf("%s\t%u\t%ld\t",divstring,chrpos,total);
368 print_line(ptr);
369 printf("\n");
370 }
371 }
372 total += subtotal;
373 ptr = nextptr;
374 if ((ptr = index(ptr,'\n')) == NULL) {
375 fprintf(stderr,"Premature end of tally from %u to %u\n",
376 Interval_low(interval),Interval_high(interval));
377 return total;
378 } else {
379 ptr++;
380 }
381 chrpos++;
382 }
383
384 *lastcoord = chrpos;
385
386 if (allocp == true) {
387 FREE(restofheader);
388 }
389
390 return total;
391 }
392
393
394 /* Need to store just the part of the query specified (e.g., 1..10) */
395 static void
print_interval_runlength(Chrpos_T * lastcoord,char * divstring,Chrpos_T coordstart,Chrpos_T coordend,int indexi,IIT_T iit,bool zeroesp)396 print_interval_runlength (Chrpos_T *lastcoord, char *divstring, Chrpos_T coordstart, Chrpos_T coordend,
397 int indexi, IIT_T iit, bool zeroesp) {
398 Interval_T interval;
399 char *label;
400 bool allocp;
401 Chrpos_T chrpos, intervalend;
402 int value;
403
404 label = IIT_label(iit,indexi,&allocp);
405 value = atoi(label);
406
407 interval = IIT_interval(iit,indexi);
408 chrpos = Interval_low(interval);
409 intervalend = Interval_high(interval);
410
411 if (zeroesp == true) {
412 while (*lastcoord < chrpos) {
413 printf("%s\t%u\t%d\n",divstring,*lastcoord,0);
414 (*lastcoord)++;
415 }
416 }
417
418 while (chrpos < coordstart) {
419 chrpos++;
420 }
421
422 while (chrpos <= intervalend && chrpos <= coordend) {
423 printf("%s\t%u\t%d\n",divstring,chrpos,value);
424 chrpos++;
425 }
426
427 *lastcoord = chrpos;
428
429 if (allocp == true) {
430 FREE(label);
431 }
432
433 return;
434 }
435
436
437 /************************************************************************
438 * Totals
439 ************************************************************************/
440
441 /* Need to store just the part of the query specified (e.g., 1..10) */
442 static void
compute_totals_tally(long int * total,int * n,Chrpos_T coordstart,Chrpos_T coordend,int indexi,IIT_T iit)443 compute_totals_tally (long int *total, int *n, Chrpos_T coordstart, Chrpos_T coordend,
444 int indexi, IIT_T iit) {
445 Interval_T interval;
446 char *annotation, *restofheader, *ptr;
447 bool allocp;
448 Chrpos_T chrpos, intervalend;
449
450 annotation = IIT_annotation(&restofheader,iit,indexi,&allocp);
451
452 interval = IIT_interval(iit,indexi);
453 chrpos = Interval_low(interval);
454 intervalend = Interval_high(interval);
455
456 ptr = annotation;
457
458 while (chrpos < coordstart) {
459 if ((ptr = index(ptr,'\n')) == NULL) {
460 fprintf(stderr,"Premature end of tally from %u to %u\n",
461 Interval_low(interval),Interval_high(interval));
462 return;
463 } else {
464 ptr++;
465 }
466 chrpos++;
467 }
468
469 while (chrpos <= intervalend && chrpos <= coordend) {
470 ptr = get_total_tally(&(*total),ptr);
471 *n += 1;
472 ptr++;
473 chrpos++;
474 }
475
476 if (allocp == true) {
477 FREE(restofheader);
478 }
479
480 return;
481 }
482
483
484 #if 0
485 /* Need to store just the part of the query specified (e.g., 1..10) */
486 static double
487 compute_logtotal_tally (long int *total, int *n, Chrpos_T coordstart, Chrpos_T coordend,
488 int indexi, IIT_T iit) {
489 double logtotal = 0.0;
490 Interval_T interval;
491 char *annotation, *restofheader, *ptr;
492 bool allocp;
493 Chrpos_T chrpos, intervalend;
494 long int count;
495
496 annotation = IIT_annotation(&restofheader,iit,indexi,&allocp);
497
498 interval = IIT_interval(iit,indexi);
499 chrpos = Interval_low(interval);
500 intervalend = Interval_high(interval);
501
502 ptr = annotation;
503
504 while (chrpos < coordstart) {
505 if ((ptr = index(ptr,'\n')) == NULL) {
506 fprintf(stderr,"Premature end of tally from %u to %u\n",
507 Interval_low(interval),Interval_high(interval));
508 return logtotal;
509 } else {
510 ptr++;
511 }
512 chrpos++;
513 }
514
515 while (chrpos <= intervalend && chrpos <= coordend) {
516 count = 0;
517 ptr = get_total_tally(&count,ptr);
518
519 logtotal += log((double) count + 1.0);
520 *total += count;
521 *n += 1;
522 ptr++;
523 chrpos++;
524 }
525
526 if (allocp == true) {
527 FREE(restofheader);
528 }
529
530 return logtotal;
531 }
532 #endif
533
534
535
536 /* coordstart used only if centerp or tallyp is true */
537 static long int
print_interval(Chrpos_T * lastcoord,long int total,char * divstring,Chrpos_T coordstart,Chrpos_T coordend,int index,IIT_T iit,int ndivs,int fieldint)538 print_interval (Chrpos_T *lastcoord, long int total,
539 char *divstring, Chrpos_T coordstart, Chrpos_T coordend,
540 int index, IIT_T iit, int ndivs, int fieldint) {
541 Interval_T interval;
542 char *label, *annotation, *restofheader, *p, *q;
543 Chrpos_T pos;
544 bool allocp;
545
546 debug(printf("Entered print_interval with coordstart %u and coordend %u\n",coordstart,coordend));
547 if (centerp == true) {
548 print_interval_centered(divstring,coordstart,index,iit,fieldint);
549 return 0;
550 } else if (tallyp == true) {
551 total += print_interval_tally(&(*lastcoord),divstring,coordstart,coordend,index,iit,zeroesp);
552 return total;
553 } else if (runlengthp == true) {
554 print_interval_runlength(&(*lastcoord),divstring,coordstart,coordend,index,iit,zeroesp);
555 return 0;
556 }
557
558 if (annotationonlyp == false) {
559 label = IIT_label(iit,index,&allocp);
560 printf(">%s ",label);
561 if (allocp == true) {
562 FREE(label);
563 }
564
565 if (ndivs > 1) {
566 if (divstring == NULL) {
567 /* For example, if interval was retrieved by label */
568 divstring = IIT_divstring_from_index(iit,index);
569 }
570 printf("%s:",divstring);
571 }
572
573 debug(printf("index is %d\n",index));
574 interval = IIT_interval(iit,index);
575 if (signedp == false) {
576 printf("%u..%u",Interval_low(interval),Interval_high(interval));
577 } else if (Interval_sign(interval) < 0) {
578 printf("%u..%u",Interval_high(interval),Interval_low(interval));
579 } else {
580 printf("%u..%u",Interval_low(interval),Interval_high(interval));
581 }
582 if (Interval_type(interval) > 0) {
583 printf(" %s",IIT_typestring(iit,Interval_type(interval)));
584 }
585 #if 0
586 /* Unnecessary because of "\n" after restofheader below */
587 if (IIT_version(iit) < 5) {
588 printf("\n");
589 }
590 #endif
591 }
592
593
594 if (fieldint < 0) {
595 annotation = IIT_annotation(&restofheader,iit,index,&allocp);
596 printf("%s\n",restofheader);
597 if (coordstart == 0 && coordend == 0) {
598 printf("%s",annotation);
599
600 } else {
601 /* printf("%.*s\n",coordend-coordstart+1,&(annotation[coordstart-1])); */
602
603 p = annotation;
604 pos = 1;
605 while (pos < coordstart) {
606 if (*p++ != '\n') {
607 pos++;
608 }
609 }
610
611 while (*p == '\n') {
612 p++;
613 }
614 q = p;
615
616 while (pos < coordend) {
617 if (*++q != '\n') {
618 pos++;
619 }
620 }
621
622 printf("%.*s\n",(int) (q-p+1),p);
623 }
624 if (allocp == true) {
625 FREE(restofheader);
626 }
627 } else {
628 annotation = IIT_annotation(&restofheader,iit,index,&allocp);
629 printf("%s\n",restofheader);
630 if (allocp == true) {
631 FREE(restofheader);
632 }
633 annotation = IIT_fieldvalue(iit,index,fieldint);
634 printf("%s\n",annotation);
635 FREE(annotation);
636 }
637
638 return 0;
639 }
640
641
642 static void
print_interval_univ(int index,Univ_IIT_T chromosome_iit)643 print_interval_univ (int index, Univ_IIT_T chromosome_iit) {
644 Univinterval_T interval;
645 char *label, *annotation, *restofheader;
646 bool allocp;
647
648 if (annotationonlyp == true) {
649 annotation = Univ_IIT_annotation(&restofheader,chromosome_iit,index,&allocp);
650 printf("%s",annotation);
651 if (allocp) {
652 FREE(restofheader);
653 }
654
655 } else {
656 label = Univ_IIT_label(chromosome_iit,index,&allocp);
657 printf(">%s ",label);
658 if (allocp == true) {
659 FREE(label);
660 }
661
662 debug(printf("index is %d\n",index));
663 interval = Univ_IIT_interval(chromosome_iit,index);
664 printf("%llu..%llu",
665 (unsigned long long) Univinterval_low(interval),
666 (unsigned long long) Univinterval_high(interval));
667 if (Univinterval_type(interval) > 0) {
668 printf(" %s",Univ_IIT_typestring(chromosome_iit,Univinterval_type(interval)));
669 }
670
671 annotation = Univ_IIT_annotation(&restofheader,chromosome_iit,index,&allocp);
672 printf("%s\n",restofheader);
673 printf("%s",annotation);
674 if (allocp) {
675 FREE(restofheader);
676 }
677
678 }
679
680 return;
681 }
682
683
684 static int *
get_matches(int * nmatches,char ** divstring,Univcoord_T * coordstart,Univcoord_T * coordend,int ** leftflanks,int * nleftflanks,int ** rightflanks,int * nrightflanks,char * query,char * typestring,IIT_T * iit,char * filename)685 get_matches (int *nmatches, char **divstring, Univcoord_T *coordstart, Univcoord_T *coordend,
686 int **leftflanks, int *nleftflanks, int **rightflanks, int *nrightflanks,
687 char *query, char *typestring, IIT_T *iit, char *filename) {
688 int *matches;
689 bool revcomp;
690 int typeint;
691
692 debug(printf("Entering get_matches with query %s.\n",query));
693
694 if (force_label_p == true || (force_coords_p == false && Parserange_query(&(*divstring),&(*coordstart),&(*coordend),&revcomp,query,filename) == false)) {
695 /* Treat query as a label */
696 *divstring = (char *) NULL;
697 *coordstart = *coordend = 0;
698
699 if (*iit == NULL) {
700 /* Read no divs */
701 if ((*iit = IIT_read(filename,/*name*/NULL,true,/*divread*/READ_NONE,/*divstring*/NULL,
702 /*add_iit_p*/true)) == NULL) {
703 if (Access_file_exists_p(filename) == false) {
704 fprintf(stderr,"Cannot read file %s\n",filename);
705 } else {
706 fprintf(stderr,"File %s appears to be an invalid IIT file\n",filename);
707 }
708 exit(9);
709
710 } else if (fieldstring != NULL) {
711 if ((fieldint = IIT_fieldint(*iit,fieldstring)) < 0) {
712 fprintf(stderr,"No field %s defined in iit file.\n",fieldstring);
713 exit(9);
714 }
715 }
716 }
717 *nleftflanks = *nrightflanks = 0;
718 matches = IIT_find(&(*nmatches),*iit,query);
719
720 } else {
721 if (*iit == NULL) {
722 if ((*iit = IIT_read(filename,/*name*/NULL,true,/*divread*/READ_ONE,*divstring,
723 /*add_iit_p*/true)) == NULL) {
724 if (Access_file_exists_p(filename) == false) {
725 fprintf(stderr,"Cannot read file %s\n",filename);
726 } else {
727 fprintf(stderr,"File %s appears to be an invalid IIT file\n",filename);
728 }
729 exit(9);
730
731 } else if (fieldstring != NULL) {
732 if ((fieldint = IIT_fieldint(*iit,fieldstring)) < 0) {
733 fprintf(stderr,"No field %s defined in iit file.\n",fieldstring);
734 exit(9);
735 }
736 }
737 }
738
739 if (typestring == NULL) {
740 /* Treat query as coordinates, without a typestring */
741 if (exactp == true) {
742 matches = IIT_get_exact_multiple(&(*nmatches),*iit,*divstring,*coordstart,*coordend,/*type*/0);
743 } else {
744 matches = IIT_get(&(*nmatches),*iit,*divstring,*coordstart,*coordend,sortp);
745 }
746 if (nflanking > 0) {
747 IIT_get_flanking(&(*leftflanks),&(*nleftflanks),&(*rightflanks),&(*nrightflanks),*iit,*divstring,
748 *coordstart,*coordend,nflanking,/*sign*/0);
749 }
750
751 } else if ((typeint = IIT_typeint(*iit,typestring)) < 0) {
752 fprintf(stderr,"No such type as %s.\n",typestring);
753 #if 0
754 /* Treat query as coordinates, without a typestring */
755 matches = IIT_get(&(*nmatches),*iit,*divstring,*coordstart,*coordend,sortp);
756 if (nflanking > 0) {
757 IIT_get_flanking(&(*leftflanks),&(*nleftflanks),&(*rightflanks),&(*nrightflanks),*iit,*divstring,
758 *coordstart,*coordend,nflanking,/*sign*/0);
759 }
760 #else
761 matches = (int *) NULL;
762 *nleftflanks = *nrightflanks = 0;
763 nmatches = 0;
764 #endif
765
766 } else {
767 /* Treat query as coordinates, with a typestring */
768 if (exactp == true) {
769 matches = IIT_get_exact_multiple(&(*nmatches),*iit,*divstring,*coordstart,*coordend,typeint);
770 } else {
771 matches = IIT_get_typed(&(*nmatches),*iit,*divstring,*coordstart,*coordend,typeint,sortp);
772 }
773 if (nflanking > 0) {
774 debug(printf("Running IIT_get_flanking_typed\n"));
775 IIT_get_flanking_typed(&(*leftflanks),&(*nleftflanks),&(*rightflanks),&(*nrightflanks),*iit,*divstring,
776 *coordstart,*coordend,nflanking,typeint,/*sign*/0);
777 }
778 }
779
780 }
781
782 return matches;
783 }
784
785 static int *
get_matches_univ(int * nmatches,char ** divstring,Univcoord_T * coordstart,Univcoord_T * coordend,char * query,char * typestring,Univ_IIT_T * chromosome_iit,char * filename)786 get_matches_univ (int *nmatches, char **divstring, Univcoord_T *coordstart, Univcoord_T *coordend,
787 char *query, char *typestring, Univ_IIT_T *chromosome_iit, char *filename) {
788 int *matches;
789 bool revcomp;
790 int typeint;
791
792 debug(printf("Entering get_matches_univ with query %s.\n",query));
793
794 if (force_label_p == true || (force_coords_p == false && Parserange_query(&(*divstring),&(*coordstart),&(*coordend),&revcomp,query,filename) == false)) {
795 /* Treat query as a label */
796 *divstring = (char *) NULL;
797 *coordstart = *coordend = 0;
798
799 if (*chromosome_iit == NULL) {
800 /* Read no divs */
801 if ((*chromosome_iit = Univ_IIT_read(filename,/*readonlyp*/true,/*add_iit_p*/true)) == NULL) {
802 if (Access_file_exists_p(filename) == false) {
803 fprintf(stderr,"Cannot read file %s\n",filename);
804 } else {
805 fprintf(stderr,"File %s appears to be an invalid IIT file\n",filename);
806 }
807 exit(9);
808 }
809 }
810 matches = Univ_IIT_find(&(*nmatches),*chromosome_iit,query);
811
812 } else {
813 if (*chromosome_iit == NULL) {
814 if ((*chromosome_iit = Univ_IIT_read(filename,/*readonlyp*/true,/*add_iit_p*/true)) == NULL) {
815 if (Access_file_exists_p(filename) == false) {
816 fprintf(stderr,"Cannot read file %s\n",filename);
817 } else {
818 fprintf(stderr,"File %s appears to be an invalid IIT file\n",filename);
819 }
820 exit(9);
821 }
822 }
823
824 if (typestring == NULL) {
825 /* Treat query as coordinates, without a typestring */
826 matches = Univ_IIT_get(&(*nmatches),*chromosome_iit,*coordstart,*coordend);
827 if (nflanking > 0) {
828 fprintf(stderr,"Flanking not supported on chromosome IIT\n");
829 }
830
831 } else if ((typeint = Univ_IIT_typeint(*chromosome_iit,typestring)) < 0) {
832 fprintf(stderr,"No such type as %s.\n",typestring);
833 #if 0
834 /* Treat query as coordinates, without a typestring */
835 matches = Univ_IIT_get(&(*nmatches),*chromosome_iit,*coordstart,*coordend);
836 if (nflanking > 0) {
837 fprintf(stderr,"Flanking not supported on chromosome IIT\n");
838 }
839 #else
840 matches = (int *) NULL;
841 *nmatches = 0;
842 #endif
843
844 } else {
845 /* Treat query as coordinates, with a typestring */
846 fprintf(stderr,"Ignoring types on chromosome IIT\n");
847 matches = Univ_IIT_get(&(*nmatches),*chromosome_iit,*coordstart,*coordend);
848 }
849
850 }
851
852 return matches;
853 }
854
855
856 static int *
get_matches_multiple_typed(int * nmatches,char ** divstring,Univcoord_T * coordstart,Univcoord_T * coordend,int ** leftflanks,int * nleftflanks,int ** rightflanks,int * nrightflanks,char * query,int * types,int ntypes,IIT_T * iit,char * filename)857 get_matches_multiple_typed (int *nmatches, char **divstring, Univcoord_T *coordstart, Univcoord_T *coordend,
858 int **leftflanks, int *nleftflanks, int **rightflanks, int *nrightflanks,
859 char *query, int *types, int ntypes, IIT_T *iit, char *filename) {
860 int *matches;
861 bool revcomp;
862
863 if (force_label_p == true || (force_coords_p == false && Parserange_query(&(*divstring),&(*coordstart),&(*coordend),&revcomp,query,filename) == false)) {
864 /* Not expecting a label */
865 abort();
866 }
867
868 if ((*iit = IIT_read(filename,/*name*/NULL,true,/*divread*/READ_ONE,*divstring,/*add_iit_p*/true)) == NULL) {
869 if (Access_file_exists_p(filename) == false) {
870 fprintf(stderr,"Cannot read file %s\n",filename);
871 } else {
872 fprintf(stderr,"File %s appears to be an invalid IIT file\n",filename);
873 }
874 exit(9);
875
876 } else if (fieldstring != NULL) {
877 if ((fieldint = IIT_fieldint(*iit,fieldstring)) < 0) {
878 fprintf(stderr,"No field %s defined in iit file.\n",fieldstring);
879 exit(9);
880 }
881 }
882
883 matches = IIT_get_multiple_typed(&(*nmatches),*iit,*divstring,*coordstart,*coordend,types,ntypes,sortp);
884 if (nflanking > 0) {
885 IIT_get_flanking_multiple_typed(&(*leftflanks),&(*nleftflanks),&(*rightflanks),&(*nrightflanks),*iit,*divstring,
886 *coordstart,*coordend,nflanking,types,ntypes);
887 }
888
889 return matches;
890 }
891
892
893 static int
int_cmp(const void * x,const void * y)894 int_cmp (const void *x, const void *y) {
895 int a = * (int *) x;
896 int b = * (int *) y;
897
898 if (a < b) {
899 return -1;
900 } else if (b < a) {
901 return +1;
902 } else {
903 return 0;
904 }
905 }
906
907
908 static int *
match_intersection(int * nmatches,int * matches1,int nmatches1,int * matches2,int nmatches2)909 match_intersection (int *nmatches, int *matches1, int nmatches1, int *matches2, int nmatches2) {
910 int *intersection;
911 int i, j, k;
912
913 qsort(matches1,nmatches1,sizeof(int),int_cmp);
914 qsort(matches2,nmatches2,sizeof(int),int_cmp);
915
916 i = j = 0;
917 k = 0;
918 while (i < nmatches1 && j < nmatches2) {
919 if (matches1[i] < matches2[j]) {
920 i++;
921 } else if (matches1[i] > matches2[j]) {
922 j++;
923 } else {
924 i++;
925 j++;
926 k++;
927 }
928 }
929
930 if ((*nmatches = k) == 0) {
931 FREE(matches1);
932 return (int *) NULL;
933
934 } else {
935 intersection = (int *) CALLOC(*nmatches,sizeof(int));
936
937 i = j = 0;
938 k = 0;
939 while (i < nmatches1 && j < nmatches2) {
940 if (matches1[i] < matches2[j]) {
941 i++;
942 } else if (matches1[i] > matches2[j]) {
943 j++;
944 } else {
945 i++;
946 j++;
947 intersection[k++] = matches1[i];
948 }
949 }
950
951 FREE(matches1);
952 return intersection;
953 }
954 }
955
956
957 int
main(int argc,char * argv[])958 main (int argc, char *argv[]) {
959 char *filename;
960 char *divstring = NULL, *lasttypestring, *ptr;
961 Univcoord_T univ_coordstart, univ_coordend;
962 Chrpos_T coordstart, coordend, lastcoord = 0U;
963 char *line, *nocomment, *query, *typestring;
964 int line_length;
965 int typeint, *types, c;
966 int nargs, ntypes, ndivs;
967 int *value_matches = NULL, *matches = NULL;
968 int n_value_matches = 0, nmatches = 0, i;
969 int *leftflanks, *rightflanks, nleftflanks = 0, nrightflanks = 0;
970 long int total;
971 int n;
972 IIT_T iit = NULL;
973 Univ_IIT_T chromosome_iit = NULL;
974 bool skipp, universalp;
975
976 int opt;
977 extern int optind;
978 extern char *optarg;
979 const char *long_name;
980 int long_option_index = 0;
981
982 while ((opt = getopt_long(argc,argv,"a:b:f:LCASUu:c:HRTZN",long_options,&long_option_index)) != -1) {
983 switch (opt) {
984 case 0:
985 long_name = long_options[long_option_index].name;
986 if (!strcmp(long_name,"version")) {
987 print_program_version();
988 exit(0);
989 } else if (!strcmp(long_name,"help")) {
990 print_program_usage();
991 exit(0);
992 } else if (!strcmp(long_name,"exact")) {
993 exactp = true;
994 } else {
995 /* Shouldn't reach here */
996 fprintf(stderr,"Don't recognize option %s. For usage, run 'gsnap --help'",long_name);
997 exit(9);
998 }
999 break;
1000
1001 case 'a': lowval = atof(optarg); user_lowvalue_p = true; value_matches_p = true; break;
1002 case 'b': highval = atof(optarg); user_highvalue_p = true; value_matches_p = true; break;
1003
1004 case 'f': fieldstring = optarg; break;
1005 case 'L': force_label_p = true; break;
1006 case 'C': force_coords_p = true; break;
1007 case 'A': annotationonlyp = true; break;
1008 case 'S': sortp = true; break;
1009 case 'U': signedp = false; break;
1010 case 'u': nflanking = atoi(optarg); break;
1011
1012 case 'c': centerp = true; centerlength = atoi(optarg); break;
1013 case 'H': centeruc = true; break;
1014 case 'R': runlengthp = true; break;
1015 case 'T': tallyp = true; break;
1016 case 'Z': zeroesp = true; break;
1017 case 'N': statsp = true; break;
1018
1019 case 'V': print_program_version(); exit(0);
1020 case '?': print_program_usage(); exit(0);
1021 default: exit(9);
1022 }
1023 }
1024
1025 argc -= (optind - 1);
1026 argv += (optind - 1);
1027
1028 if (argc <= 1) {
1029 fprintf(stderr,"Need to specify an iit file. Type \"iit_get --help\" for help.\n");
1030 exit(9);
1031 } else {
1032 filename = argv[1];
1033 }
1034
1035 if (value_matches_p == true) {
1036 if ((iit = IIT_read(filename,/*name*/NULL,true,/*divread*/READ_ALL,/*divstring*/NULL,
1037 /*add_iit_p*/true)) == NULL) {
1038 if (Access_file_exists_p(filename) == false) {
1039 fprintf(stderr,"Cannot read file %s\n",filename);
1040 } else {
1041 fprintf(stderr,"File %s appears to be an invalid IIT file\n",filename);
1042 }
1043 exit(9);
1044
1045 } else if (fieldstring != NULL) {
1046 if ((fieldint = IIT_fieldint(iit,fieldstring)) < 0) {
1047 fprintf(stderr,"No field %s defined in iit file.\n",fieldstring);
1048 exit(9);
1049 }
1050 }
1051
1052 if (IIT_valuep(iit) == false) {
1053 fprintf(stderr,"Error: This IIT file does not have values stored\n");
1054 exit(9);
1055 }
1056
1057 if (user_lowvalue_p == true && user_highvalue_p == true) {
1058 if (lowval > highval) {
1059 fprintf(stderr,"Cannot have lowval %f > highval %f\n",lowval,highval);
1060 exit(9);
1061 } else {
1062 value_matches = IIT_get_values_between(&n_value_matches,iit,lowval,highval);
1063 }
1064 } else if (user_lowvalue_p == true) {
1065 value_matches = IIT_get_values_above(&n_value_matches,iit,lowval);
1066
1067 } else { /* user_highvalue_p == true */
1068 value_matches = IIT_get_values_below(&n_value_matches,iit,highval);
1069 }
1070 }
1071
1072 if (0 && statsp == true && argc == 2) {
1073 /* Want total over entire IIT */
1074 if ((iit = IIT_read(filename,NULL,true,/*divread*/READ_ALL,/*divstring*/NULL,/*add_iit_p*/true)) == NULL) {
1075 if (Access_file_exists_p(filename) == false) {
1076 fprintf(stderr,"Cannot read file %s\n",filename);
1077 } else {
1078 fprintf(stderr,"File %s appears to be an invalid IIT file\n",filename);
1079 }
1080 exit(9);
1081
1082 } else if (fieldstring != NULL) {
1083 if ((fieldint = IIT_fieldint(iit,fieldstring)) < 0) {
1084 fprintf(stderr,"No field %s defined in iit file.\n",fieldstring);
1085 exit(9);
1086 }
1087 }
1088
1089 total = 0;
1090 n = 0;
1091 for (i = 0; i < IIT_total_nintervals(iit); i++) {
1092 debug(printf("index = %d\n",matches[i]));
1093 compute_totals_tally(&total,&n,/*coordstart*/0,/*coordend*/-1U,i,iit);
1094 }
1095 printf("counts:%ld non-zero-positions:%u mean-over-nonzero:%.3f\n",total,n,(double) total/(double) n);
1096
1097 IIT_free(&iit);
1098 return 0;
1099
1100 } else if ((universalp = IIT_universalp(filename,/*add_iit_p*/true)) == true) {
1101 chromosome_iit = Univ_IIT_read(filename,/*readonlyp*/true,/*add_iit_p*/true);
1102 if (argc != 3) {
1103 fprintf(stderr,"For chromosome IIT file, need to specify a query on the command line\n");
1104 exit(9);
1105 } else {
1106 /* Try as 0:<iitfile> 1:<query> */
1107 matches = get_matches_univ(&nmatches,&divstring,&univ_coordstart,&univ_coordend,
1108 argv[2],/*typestring*/NULL,&chromosome_iit,filename);
1109
1110 for (i = 0; i < nmatches; i++) {
1111 debug(printf("\nindex = %d\n",matches[i]));
1112 print_interval_univ(matches[i],chromosome_iit);
1113 }
1114 FREE(matches);
1115 }
1116
1117 if (divstring != NULL) {
1118 FREE(divstring);
1119 }
1120 Univ_IIT_free(&chromosome_iit);
1121 return 0;
1122
1123 } else if (argc == 2 && value_matches_p == true) {
1124 /* Note: Could potentially handle input from stdin, but currently just deal with value_matches */
1125
1126 ndivs = IIT_ndivs(iit);
1127 for (i = 0; i < n_value_matches; i++) {
1128 debug(printf("\nindex = %d\n",matches[i]));
1129 print_interval(&lastcoord,/*total*/0,/*divstring*/NULL,/*coordstart*/0,/*coordend*/0,
1130 value_matches[i],iit,ndivs,fieldint);
1131 }
1132
1133 FREE(value_matches);
1134 IIT_free(&iit);
1135 return 0;
1136
1137 } else if (argc == 2) {
1138 debug(printf("Running under argc 2\n"));
1139 /* Expecting input from stdin */
1140
1141 if ((iit = IIT_read(filename,NULL,true,/*divread*/READ_ALL,/*divstring*/NULL,/*add_iit_p*/true)) == NULL) {
1142 if (Access_file_exists_p(filename) == false) {
1143 fprintf(stderr,"Cannot read file %s\n",filename);
1144 } else {
1145 fprintf(stderr,"File %s appears to be an invalid IIT file\n",filename);
1146 }
1147 exit(9);
1148 } else if (fieldstring != NULL) {
1149 if ((fieldint = IIT_fieldint(iit,fieldstring)) < 0) {
1150 fprintf(stderr,"No field %s defined in iit file.\n",fieldstring);
1151 exit(9);
1152 }
1153 }
1154
1155 while ((line = Getline_wlength(&line_length,stdin)) != NULL) {
1156 #if 0
1157 if ((ptr = rindex(Buffer,'\n')) != NULL) {
1158 *ptr = '\0';
1159 }
1160 #endif
1161 nocomment = (char *) MALLOC((line_length+1)*sizeof(char));
1162 strcpy(nocomment,line);
1163
1164 if ((ptr = rindex(nocomment,'#')) != NULL) {
1165 *ptr = '\0';
1166 }
1167
1168 skipp = false;
1169
1170 divstring = (char *) NULL;
1171 query = (char *) MALLOC((line_length+1)*sizeof(char));
1172 typestring = (char *) MALLOC((line_length+1)*sizeof(char));
1173
1174 if ((nargs = sscanf(nocomment,"%s %s",query,typestring)) == 2) {
1175 debug(printf("typestring is %s\n",typestring));
1176 matches = get_matches(&nmatches,&divstring,&univ_coordstart,&univ_coordend,
1177 &leftflanks,&nleftflanks,&rightflanks,&nrightflanks,
1178 query,typestring,&iit,filename);
1179 coordstart = (Chrpos_T) univ_coordstart;
1180 coordend = (Chrpos_T) univ_coordend;
1181
1182 } else if (nargs == 1) {
1183 debug(printf("typestring is NULL\n"));
1184 matches = get_matches(&nmatches,&divstring,&univ_coordstart,&univ_coordend,
1185 &leftflanks,&nleftflanks,&rightflanks,&nrightflanks,
1186 query,/*typestring*/NULL,&iit,filename);
1187 coordstart = (Chrpos_T) univ_coordstart;
1188 coordend = (Chrpos_T) univ_coordend;
1189
1190 } else {
1191 fprintf(stderr,"Can't parse line %s. Ignoring.\n",nocomment);
1192 skipp = true;
1193 }
1194 FREE(divstring);
1195
1196 total = 0;
1197 if (skipp == false) {
1198 fprintf(stdout,"# Query: %s\n",line);
1199 ndivs = IIT_ndivs(iit);
1200 if (nflanking > 0) {
1201 for (i = nleftflanks-1; i >= 0; i--) {
1202 debug(printf("\nleft index = %d\n",leftflanks[i]));
1203 print_interval(&lastcoord,/*total*/0,divstring,coordstart,coordend,leftflanks[i],iit,ndivs,fieldint);
1204 }
1205 printf("====================\n");
1206 FREE(leftflanks);
1207 }
1208
1209 lastcoord = coordstart;
1210 for (i = 0; i < nmatches; i++) {
1211 debug(printf("\nindex = %d\n",matches[i]));
1212 total = print_interval(&lastcoord,total,divstring,coordstart,coordend,matches[i],iit,ndivs,fieldint);
1213 }
1214
1215 if (nflanking > 0) {
1216 printf("====================\n");
1217 for (i = 0; i < nrightflanks; i++) {
1218 debug(printf("\nright index = %d\n",rightflanks[i]));
1219 print_interval(&lastcoord,/*total*/0,divstring,coordstart,coordend,rightflanks[i],iit,ndivs,fieldint);
1220 }
1221 FREE(rightflanks);
1222 }
1223
1224 if (zeroesp == true) {
1225 while (lastcoord <= coordend) {
1226 printf("%s\t%u\t%d\n",divstring,lastcoord,0);
1227 lastcoord++;
1228 }
1229 }
1230 }
1231
1232 FREE(divstring);
1233 FREE(matches);
1234 /* printf("%ld\n",total); -- Not sure why this was here */
1235 fprintf(stdout,"# End\n");
1236 fflush(stdout);
1237
1238 FREE(typestring);
1239 FREE(query);
1240 FREE(nocomment);
1241 FREE(line);
1242 }
1243
1244 IIT_free(&iit);
1245 return 0;
1246
1247 } else {
1248 /* Get coordinates/type from command line */
1249 divstring = (char *) NULL;
1250 if (argc == 3) {
1251 /* Try as 0:<iitfile> 1:<query> */
1252 matches = get_matches(&nmatches,&divstring,&univ_coordstart,&univ_coordend,
1253 &leftflanks,&nleftflanks,&rightflanks,&nrightflanks,
1254 argv[2],/*typestring*/NULL,&iit,filename);
1255 coordstart = (Chrpos_T) univ_coordstart;
1256 coordend = (Chrpos_T) univ_coordend;
1257
1258 } else if (argc == 4) {
1259 /* Try as 0:<iitfile> 1:<query> 2:<type> */
1260 divstring = (char *) NULL;
1261 debug(printf("Running under argc 4\n"));
1262 matches = get_matches(&nmatches,&divstring,&univ_coordstart,&univ_coordend,
1263 &leftflanks,&nleftflanks,&rightflanks,&nrightflanks,
1264 argv[2],argv[3],&iit,filename);
1265 coordstart = (Chrpos_T) univ_coordstart;
1266 coordend = (Chrpos_T) univ_coordend;
1267
1268 } else {
1269 types = (int *) CALLOC(argc-3,sizeof(int));
1270 for (c = 3, ntypes = 0; c < argc; c++) {
1271 if ((typeint = IIT_typeint(iit,argv[c])) < 0) {
1272 fprintf(stderr,"No such type as %s. Ignoring the type.\n",argv[c]);
1273 } else {
1274 types[ntypes++] = typeint;
1275 lasttypestring = argv[c];
1276 }
1277 }
1278 if (ntypes == 0) {
1279 matches = get_matches(&nmatches,&divstring,&univ_coordstart,&univ_coordend,
1280 &leftflanks,&nleftflanks,&rightflanks,&nrightflanks,
1281 argv[2],/*typestring*/NULL,&iit,filename);
1282 coordstart = (Chrpos_T) univ_coordstart;
1283 coordend = (Chrpos_T) univ_coordend;
1284
1285 } else if (ntypes == 1) {
1286 matches = get_matches(&nmatches,&divstring,&univ_coordstart,&univ_coordend,
1287 &leftflanks,&nleftflanks,&rightflanks,&nrightflanks,
1288 argv[2],lasttypestring,&iit,filename);
1289 coordstart = (Chrpos_T) univ_coordstart;
1290 coordend = (Chrpos_T) univ_coordend;
1291
1292 } else {
1293 matches = get_matches_multiple_typed(&nmatches,&divstring,&univ_coordstart,&univ_coordend,
1294 &leftflanks,&nleftflanks,&rightflanks,&nrightflanks,
1295 argv[2],types,ntypes,&iit,filename);
1296 coordstart = (Chrpos_T) univ_coordstart;
1297 coordend = (Chrpos_T) univ_coordend;
1298 }
1299 }
1300
1301 if (value_matches_p == true) {
1302 matches = match_intersection(&nmatches,/*matches1*/matches,/*nmatches1*/nmatches,
1303 /*matches2*/value_matches,/*nmatches2*/n_value_matches);
1304 FREE(value_matches);
1305 }
1306
1307 #if 0
1308 if (centerp == true) {
1309 print_spaces(centerlength);
1310 printf("*");
1311 print_spaces(centerlength-1);
1312 printf("\n");
1313 }
1314 #endif
1315
1316 if (statsp == true) {
1317 total = 0;
1318 n = 0;
1319 for (i = 0; i < nmatches; i++) {
1320 debug(printf("index = %d\n",matches[i]));
1321 compute_totals_tally(&total,&n,coordstart,coordend,matches[i],iit);
1322 }
1323 n = coordend - coordstart + 1;
1324 printf("counts:%ld width:%u mean:%.3f\n",total,n,(double)total/(double) n);
1325
1326 #if 0
1327 } else if (geomeanp == true) {
1328 logtotal = 0.0;
1329 total = 0;
1330 n = 0;
1331 for (i = 0; i < nmatches; i++) {
1332 debug(printf("index = %d\n",matches[i]));
1333 logtotal = compute_logtotal_tally(&total,&n,coordstart,coordend,matches[i],iit);
1334 }
1335 printf("geomean:%f totalcounts:%ld posrange:%d\n",
1336 exp(logtotal/(double) (coordend - coordstart + 1)) - 1.0,total,n);
1337 #endif
1338
1339 } else {
1340 ndivs = IIT_ndivs(iit);
1341 if (nflanking > 0) {
1342 for (i = nleftflanks-1; i >= 0; i--) {
1343 debug(printf("\nleft index = %d\n",leftflanks[i]));
1344 print_interval(&lastcoord,/*total*/0,divstring,coordstart,coordend,leftflanks[i],iit,ndivs,fieldint);
1345 }
1346 printf("====================\n");
1347 FREE(leftflanks);
1348 }
1349
1350 lastcoord = coordstart;
1351 for (i = 0; i < nmatches; i++) {
1352 debug(printf("\nindex = %d\n",matches[i]));
1353 print_interval(&lastcoord,/*total*/0,divstring,coordstart,coordend,matches[i],iit,ndivs,fieldint);
1354 }
1355
1356 if (nflanking > 0) {
1357 printf("====================\n");
1358 for (i = 0; i < nrightflanks; i++) {
1359 debug(printf("\nright index = %d\n",rightflanks[i]));
1360 print_interval(&lastcoord,/*total*/0,divstring,coordstart,coordend,rightflanks[i],iit,ndivs,fieldint);
1361 }
1362 FREE(rightflanks);
1363 }
1364 }
1365
1366 FREE(divstring);
1367 FREE(matches);
1368 IIT_free(&iit);
1369 return 0;
1370 }
1371 }
1372
1373