[Bio] / Kmers2 / kguts.cc Repository:
ViewVC logotype

Annotation of /Kmers2/kguts.cc

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.1 - (view) (download) (as text)

1 : olson 1.1 #include "kguts.h"
2 :    
3 :     static const char genetic_code[64] = {
4 :     'K','N','K','N','T','T','T','T','R','S','R','S','I','I','M','I',
5 :     'Q','H','Q','H','P','P','P','P','R','R','R','R','L','L','L','L',
6 :     'E','D','E','D','A','A','A','A','G','G','G','G','V','V','V','V',
7 :     '*','Y','*','Y','S','S','S','S','*','C','W','C','L','F','L','F'
8 :     };
9 :    
10 :     static const char prot_alpha[20] = {
11 :     'A','C','D','E','F','G','H','I','K','L','M','N','P','Q','R','S','T','V','W','Y'
12 :     };
13 :    
14 :    
15 :    
16 :     KmerGuts::KmerGuts()
17 :     {
18 :     tot_lookups = 0;
19 :     retry = 0;
20 :    
21 :     debug = 0;
22 :     aa = 0;
23 :     hits_only = 0;
24 :     size_hash = 1400303159;
25 :     write_mem_map = 0;
26 :     data_dir = 0;
27 :    
28 :     num_hits = 0;
29 :     num_oI = 0;
30 :    
31 :     order_constraint = 0;
32 :     min_hits = 5;
33 :     min_weighted_hits = 0;
34 :     max_gap = 200;
35 :    
36 :     pIseq = (unsigned char *) malloc(MAX_SEQ_LEN / 3);
37 :     cdata = (char *) malloc(MAX_SEQ_LEN);
38 :     data = (char *) malloc(MAX_SEQ_LEN);
39 :     pseq = (char *) malloc(MAX_SEQ_LEN / 3);
40 :     }
41 :    
42 :    
43 :     /* =========================== end of reduction global variables ================= */
44 :    
45 :     unsigned char KmerGuts::to_amino_acid_off(char c) {
46 :     switch (c)
47 :     {
48 :     case 'A':
49 :     return 0;
50 :    
51 :     case 'C':
52 :     return 1;
53 :    
54 :     case 'D':
55 :     return 2;
56 :    
57 :     case 'E':
58 :     return 3;
59 :    
60 :     case 'F':
61 :     return 4;
62 :    
63 :     case 'G':
64 :     return 5;
65 :    
66 :     case 'H':
67 :     return 6;
68 :    
69 :     case 'I':
70 :     return 7;
71 :    
72 :     case 'K':
73 :     return 8;
74 :    
75 :     case 'L':
76 :     return 9;
77 :    
78 :     case 'M':
79 :     return 10;
80 :    
81 :     case 'N':
82 :     return 11;
83 :    
84 :     case 'P':
85 :     return 12;
86 :    
87 :     case 'Q':
88 :     return 13;
89 :    
90 :     case 'R':
91 :     return 14;
92 :    
93 :     case 'S':
94 :     return 15;
95 :    
96 :     case 'T':
97 :     return 16;
98 :    
99 :     case 'V':
100 :     return 17;
101 :    
102 :     case 'W':
103 :     return 18;
104 :    
105 :     case 'Y':
106 :     return 19;
107 :    
108 :     default:
109 :     return 20;
110 :     }
111 :     }
112 :    
113 :     char KmerGuts::comp(char c)
114 :     {
115 :     switch (c)
116 :     {
117 :     case 'a':
118 :     return 't';
119 :     case 'A':
120 :     return 'T';
121 :    
122 :     case 'c':
123 :     return 'g';
124 :     case 'C':
125 :     return 'G';
126 :    
127 :     case 'g':
128 :     return 'c';
129 :     case 'G':
130 :     return 'C';
131 :    
132 :     case 't':
133 :     case 'u':
134 :     return 'a';
135 :     case 'T':
136 :     case 'U':
137 :     return 'A';
138 :    
139 :     case 'm':
140 :     return 'k';
141 :     case 'M':
142 :     return 'K';
143 :    
144 :     case 'r':
145 :     return 'y';
146 :     case 'R':
147 :     return 'Y';
148 :    
149 :     case 'w':
150 :     return 'w';
151 :     case 'W':
152 :     return 'W';
153 :    
154 :     case 's':
155 :     return 'S';
156 :     case 'S':
157 :     return 'S';
158 :    
159 :     case 'y':
160 :     return 'r';
161 :     case 'Y':
162 :     return 'R';
163 :    
164 :     case 'k':
165 :     return 'm';
166 :     case 'K':
167 :     return 'M';
168 :    
169 :     case 'b':
170 :     return 'v';
171 :     case 'B':
172 :     return 'V';
173 :    
174 :     case 'd':
175 :     return 'h';
176 :     case 'D':
177 :     return 'H';
178 :    
179 :     case 'h':
180 :     return 'd';
181 :     case 'H':
182 :     return 'D';
183 :    
184 :     case 'v':
185 :     return 'b';
186 :     case 'V':
187 :     return 'B';
188 :    
189 :     case 'n':
190 :     return 'n';
191 :     case 'N':
192 :     return 'N';
193 :    
194 :     default:
195 :     return c;
196 :     }
197 :     }
198 :    
199 :     void KmerGuts::rev_comp(char *data,char *cdata) {
200 :    
201 :     int n = strlen(data);
202 :     char *p = data + (n-1);
203 :     char *pc = cdata;
204 :     while (n--) {
205 :     *(pc++) = compl(*(p--));
206 :     }
207 :     *pc = 0;
208 :     }
209 :    
210 :     unsigned long long KmerGuts::encoded_kmer(unsigned char *p) {
211 :     unsigned long long encodedK = *p;
212 :     int i;
213 :     for (i=1; (i <= K-1); i++) {
214 :     encodedK = (encodedK * 20) + *(p+i);
215 :     }
216 :    
217 :     if (encodedK > MAX_ENCODED) {
218 :     fprintf(stderr,"bad encoding - input must have included invalid characters\n");
219 :     for (i=0; (i < K); i++) {
220 :     fprintf(stderr,"%d ",*(p+i));
221 :     }
222 :     fprintf(stderr,"\n");
223 :     exit(2);
224 :     }
225 :     return encodedK;
226 :     }
227 :    
228 :     unsigned long long KmerGuts::encoded_aa_kmer(char *p)
229 :     {
230 :     unsigned char aa_off[K];
231 :     int j;
232 :     for (j=0; (j < K); j++) {
233 :     int prot_c = *(p+j);
234 :     aa_off[j] = to_amino_acid_off(prot_c);
235 :     }
236 :     return encoded_kmer(aa_off);
237 :     }
238 :    
239 :     void KmerGuts::decoded_kmer(unsigned long long encodedK,char *decoded) {
240 :    
241 :     int i;
242 :     *(decoded+K) = '\0';
243 :     unsigned long long x = encodedK;
244 :    
245 :     for (i=K-1; (i >= 0); i--) {
246 :     *(decoded+i) = prot_alpha[x % 20];
247 :     x = x / 20;
248 :     }
249 :     }
250 :    
251 :    
252 :     int KmerGuts::dna_char(char c)
253 :     {
254 :     switch (c)
255 :     {
256 :     case 'a':
257 :     case 'A':
258 :     return 0;
259 :    
260 :     case 'c':
261 :     case 'C':
262 :     return 1;
263 :    
264 :     case 'g':
265 :     case 'G':
266 :     return 2;
267 :    
268 :     case 't':
269 :     case 'u':
270 :     case 'T':
271 :     case 'U':
272 :     return 3;
273 :    
274 :     default:
275 :     return 4;;
276 :     }
277 :     }
278 :    
279 :     void KmerGuts::translate(char *seq,int off,char *pseq, unsigned char *pIseq) {
280 :    
281 :     int i;
282 :     int max = strlen(seq) - 3;
283 :     char *p = pseq;
284 :     unsigned char *pI = pIseq;
285 :     for (i=off; (i <= max); ) {
286 :     int c1 = dna_char(seq[i++]);
287 :     int c2 = dna_char(seq[i++]);
288 :     int c3 = dna_char(seq[i++]);
289 :     if ((c1 < 4) && (c2 < 4) && (c3 < 4)) {
290 :     int I = (c1 * 16) + (c2 * 4) + c3;
291 :     char prot_c = genetic_code[I];
292 :     *(p++) = prot_c;
293 :     *(pI++) = to_amino_acid_off(prot_c);
294 :     }
295 :     else {
296 :     *(p++) = 'x';
297 :     *(pI++) = 20;
298 :     }
299 :     }
300 :     *p = 0;
301 :     *pI = 21;
302 :     if (debug >= 3) {
303 :     fprintf(stderr,"len-seq=%d max=%d p=%ld\n",(int) strlen(seq),max,p-pseq);
304 :     }
305 :     }
306 :    
307 :     #define MAX_FUNC_OI_INDEX 1000000
308 :     #define MAX_FUNC_OI_VALS 100000000
309 :    
310 :     char **KmerGuts::load_indexed_ar(char *filename,int *sz) {
311 :     char **index_ar = (char **) malloc(MAX_FUNC_OI_INDEX * sizeof(char *));
312 :     char *vals = (char *)malloc(MAX_FUNC_OI_VALS);
313 :     char *p = vals;
314 :     FILE *ifp = fopen(filename,"r");
315 :     if (ifp == NULL) {
316 :     fprintf(stderr,"could not open %s\n",filename);
317 :     exit(1);
318 :     }
319 :    
320 :     *sz = 0;
321 :     int j;
322 :     while ((fscanf(ifp,"%d\t",&j) == 1) && fgets(p,1000,ifp)) {
323 :     if (*sz != j) {
324 :     fprintf(stderr,"Your index must be dense and in order (see line %ld, should be %d)\n",p-vals,*sz);
325 :     exit(1);
326 :     }
327 :     /* fprintf(stderr,"%d is %s\n",*sz,p); */
328 :     index_ar[*sz] = p; /* the fgets leaves the \n at the end of each line */
329 :     p += strlen(index_ar[*sz]) -1;
330 :     *(p++) = '\0';
331 :     if ((*sz >= MAX_FUNC_OI_INDEX) || ((p-vals) > (MAX_FUNC_OI_VALS - 1000))) {
332 :     fprintf(stderr,"Your function or oI index arrays are too small; bump MAX_FUNC_OI_INDEX and MAX_FUNC_OI_VALS\n");
333 :     exit(1);
334 :     }
335 :    
336 :     *sz += 1;
337 :     }
338 :     return index_ar;
339 :     }
340 :    
341 :     char **KmerGuts::load_functions(char *file) {
342 :     int sz;
343 :     return load_indexed_ar(file,&sz);
344 :     }
345 :    
346 :     char **KmerGuts::load_otus(char *file) {
347 :     int sz;
348 :     return load_indexed_ar(file,&sz);
349 :     }
350 :    
351 :     long long KmerGuts::find_empty_hash_entry(sig_kmer_t sig_kmers[],unsigned long long encodedK) {
352 :     long long hash_entry = encodedK % size_hash;
353 :     while (sig_kmers[hash_entry].which_kmer <= MAX_ENCODED)
354 :     hash_entry = (hash_entry+1)%size_hash;
355 :     return hash_entry;
356 :     }
357 :    
358 :     long long KmerGuts::lookup_hash_entry(sig_kmer_t sig_kmers[],unsigned long long encodedK) {
359 :     long long hash_entry = encodedK % size_hash;
360 :     if (debug >= 2)
361 :     tot_lookups++;
362 :     while ((sig_kmers[hash_entry].which_kmer <= MAX_ENCODED) && (sig_kmers[hash_entry].which_kmer != encodedK)) {
363 :     if (debug >= 2)
364 :     retry++;
365 :     hash_entry++;
366 :     if (hash_entry == size_hash)
367 :     hash_entry = 0;
368 :     }
369 :     if (sig_kmers[hash_entry].which_kmer > MAX_ENCODED) {
370 :     return -1;
371 :     }
372 :     else {
373 :     return hash_entry;
374 :     }
375 :     }
376 :    
377 :     KmerGuts::kmer_memory_image_t *KmerGuts::load_raw_kmers(char *file,unsigned long long num_entries, unsigned long long *alloc_sz) {
378 :     /*
379 :     * Allocate enough memory to hold the kmer_memory_image_t header plus the hash table itself.
380 :     */
381 :     *alloc_sz = sizeof(kmer_memory_image_t) + (sizeof(sig_kmer_t) * num_entries);
382 :    
383 :     kmer_memory_image_t *image = (kmer_memory_image_t *) malloc(*alloc_sz);
384 :    
385 :     /*
386 :     * Initialize our table pointer to the first byte after the header.
387 :     */
388 :     image->num_sigs = num_entries;
389 :     image->entry_size = sizeof(sig_kmer_t);
390 :     image->version = (long long) VERSION;
391 :    
392 :     sig_kmer_t *sig_kmers = (sig_kmer_t *) (image + 1);
393 :    
394 :     FILE *ifp = fopen(file,"r");
395 :     if (ifp == NULL) {
396 :     fprintf(stderr,"could not open %s",file);
397 :     exit(1);
398 :     }
399 :    
400 :     long long i;
401 :     for (i=0; (i < size_hash); i++)
402 :     sig_kmers[i].which_kmer = MAX_ENCODED + 1;
403 :    
404 :     char kmer_string[K+1];
405 :     int end_off;
406 :     int fI;
407 :     float f_wt;
408 :     int oI;
409 :     long long loaded = 0;
410 :     while (fscanf(ifp,"%s\t%d\t%d\t%f\t%d",
411 :     kmer_string,&end_off,&fI,&f_wt,&oI) >= 4) {
412 :     unsigned long long encodedK = encoded_aa_kmer(kmer_string);
413 :     long long hash_entry = find_empty_hash_entry(sig_kmers,encodedK);
414 :     loaded++;
415 :     if (loaded >= (size_hash / 2)) {
416 :     fprintf(stderr,"Your Kmer hash is half-full; use -s (and -w) to bump it\n");
417 :     exit(1);
418 :     }
419 :     sig_kmers[hash_entry].which_kmer = encodedK;
420 :     sig_kmers[hash_entry].avg_from_end = end_off;
421 :     sig_kmers[hash_entry].function_index = fI;
422 :     sig_kmers[hash_entry].otu_index = oI;
423 :     sig_kmers[hash_entry].function_wt = f_wt;
424 :     }
425 :     if (debug >= 2)
426 :     fprintf(stderr,"loaded %lld kmers\n",loaded);
427 :    
428 :     return image;
429 :     }
430 :    
431 :     KmerGuts::kmer_handle_t *KmerGuts::init_kmers(char *dataD) {
432 :     kmer_handle_t *handle = (kmer_handle_t *) malloc(sizeof(kmer_handle_t));
433 :    
434 :     kmer_memory_image_t *image;
435 :    
436 :     char file[300];
437 :     strcpy(file,dataD);
438 :     strcat(file,"/function.index");
439 :     handle->function_array = load_functions(file);
440 :    
441 :     strcpy(file,dataD);
442 :     strcat(file,"/otu.index");
443 :     handle->otu_array = load_otus(file);
444 :    
445 :     char fileM[300];
446 :     strcpy(fileM,dataD);
447 :     strcat(fileM,"/kmer.table.mem_map");
448 :    
449 :     if (write_mem_map) {
450 :     unsigned long long sz, table_size;
451 :     strcpy(file,dataD);
452 :     strcat(file,"/final.kmers");
453 :    
454 :     unsigned long long image_size;
455 :    
456 :     image = load_raw_kmers(file, size_hash, &image_size);
457 :    
458 :     handle->kmer_table = (sig_kmer_t *) (image + 1);
459 :     handle->num_sigs = image->num_sigs;
460 :    
461 :     FILE *fp = fopen(fileM,"w");
462 :     if (fp == NULL) {
463 :     fprintf(stderr,"could not open %s for writing: %s ",fileM, strerror(errno));
464 :     exit(1);
465 :     }
466 :     fwrite(image, image_size, 1, fp);
467 :     fclose(fp);
468 :    
469 :     strcpy(fileM,dataD);
470 :     strcat(fileM,"/size_hash.and.table_size");
471 :     fp = fopen(fileM,"w");
472 :     fprintf(fp,"%lld\t%lld\n",sz,table_size);
473 :     fclose(fp);
474 :     }
475 :     else {
476 :     int fd;
477 :     if ((fd = open(fileM, O_RDONLY)) == -1) {
478 :     perror("open");
479 :     exit(1);
480 :     }
481 :    
482 :     /*
483 :     * Set up for creating memory image from file. Start by determining file size
484 :     * on disk with a stat() call.
485 :     */
486 :     struct stat sbuf;
487 :     if (stat(fileM, &sbuf) == -1) {
488 :     fprintf(stderr, "stat %s failed: %s\n", fileM, strerror(errno));
489 :     exit(1);
490 :     }
491 :     unsigned long long file_size = sbuf.st_size;
492 :    
493 :     /*
494 :     * Memory map.
495 :     */
496 :     int flags = MAP_SHARED;
497 :     #ifdef MAP_POPULATE
498 :     flags |= MAP_POPULATE;
499 :     #endif
500 :    
501 :     image = (kmer_memory_image_t *) mmap((caddr_t)0, file_size, PROT_READ, flags, fd, 0);
502 :    
503 :     if (image == (kmer_memory_image_t *)(-1)) {
504 :     fprintf(stderr, "mmap of kmer_table %s failed: %s\n", fileM, strerror(errno));
505 :     exit(1);
506 :     }
507 :    
508 :     /*
509 :     * Our image is mapped. Validate against the current version of this code.
510 :     */
511 :     if (image->version != (long long) VERSION) {
512 :     fprintf(stderr, "Version mismatch for file %s: file has %lld code has %lld\n",
513 :     fileM, image->version, (long long) VERSION);
514 :     exit(1);
515 :     }
516 :    
517 :     if (image->entry_size != (unsigned long long) sizeof(sig_kmer_t)) {
518 :     fprintf(stderr, "Version mismatch for file %s: file has entry size %lld code has %lld\n",
519 :     fileM, image->entry_size, (unsigned long long) sizeof(sig_kmer_t));
520 :     exit(1);
521 :     }
522 :    
523 :     size_hash = image->num_sigs;
524 :     handle->num_sigs = size_hash;
525 :     handle->kmer_table = (sig_kmer_t *) (image + 1);
526 :    
527 :     /* Validate overall file size vs the entry size and number of entries */
528 :     if (file_size != ((sizeof(sig_kmer_t) * image->num_sigs) + sizeof(kmer_memory_image_t))) {
529 :     fprintf(stderr, "Version mismatch for file %s: file size does not match\n", fileM);
530 :     exit(1);
531 :     }
532 :    
533 :     fprintf(stderr, "Set size_hash=%lld from file size %lld\n", size_hash, file_size);
534 :    
535 :     }
536 :     return handle;
537 :     }
538 :    
539 :     void KmerGuts::advance_past_ambig(unsigned char **p,unsigned char *bound) {
540 :    
541 :     if (K == 5) {
542 :     while (((*p) < bound) &&
543 :     ((*(*p) == 20) ||
544 :     (*((*p)+1) == 20) ||
545 :     (*((*p)+2) == 20) ||
546 :     (*((*p)+3) == 20) ||
547 :     (*((*p)+4) == 20) )) {
548 :     (*p)++;
549 :     }
550 :     }
551 :     else { /* ##### ASSUMING K == 8 #### */
552 :     int bad = 1;
553 :     while ((*p < bound) && (bad == 1)) {
554 :     bad = 0;
555 :     if (*((*p)+7) == 20) {
556 :     bad = 1;
557 :     (*p) += 8;
558 :     }
559 :     else if (*((*p)+6) == 20) {
560 :     bad = 1;
561 :     (*p) += 7;
562 :     }
563 :     else if (*((*p)+5) == 20) {
564 :     bad = 1;
565 :     (*p) += 6;
566 :     }
567 :     else if (*((*p)+4) == 20) {
568 :     bad = 1;
569 :     (*p) += 5;
570 :     }
571 :     else if (*((*p)+3) == 20) {
572 :     bad = 1;
573 :     (*p) += 4;
574 :     }
575 :     else if (*((*p)+2) == 20) {
576 :     bad = 1;
577 :     (*p) += 3;
578 :     }
579 :     else if (*((*p)+1) == 20) {
580 :     bad = 1;
581 :     (*p) += 2;
582 :     }
583 :     else if (*((*p)+0) == 20) {
584 :     bad = 1;
585 :     (*p) += 1;
586 :     }
587 :     }
588 :     }
589 :     }
590 :    
591 :     void KmerGuts::display_hits(FILE *fh) {
592 :     fprintf(fh, "hits: ");
593 :     int i;
594 :     for (i=0; (i < num_hits); i++) {
595 :     fprintf(fh, "%d/%f/%d ", hits[i].from0_in_prot,hits[i].function_wt,hits[i].fI);
596 :     }
597 :     fprintf(fh, "\n");
598 :     }
599 :    
600 :    
601 :     void KmerGuts::process_set_of_hits(kmer_handle_t *kmersH, FILE *fh) {
602 :     int fI_count = 0;
603 :     float weighted_hits = 0;
604 :     int last_hit=0;
605 :     int i=0;
606 :     while (i < num_hits) {
607 :     if (hits[i].fI == current_fI) {
608 :     last_hit = i;
609 :     fI_count++;
610 :     weighted_hits += hits[i].function_wt;
611 :     }
612 :     i++;
613 :     } if ((fI_count >= min_hits) && (weighted_hits >= min_weighted_hits)) {
614 :     if (!hits_only)
615 :     fprintf(fh, "CALL\t%d\t%d\t%d\t%d\t%s\t%f\n",hits[0].from0_in_prot,
616 :     hits[last_hit].from0_in_prot+(K-1),
617 :     fI_count,
618 :     current_fI,
619 :     kmersH->function_array[current_fI],
620 :     weighted_hits);
621 :    
622 :     if (debug > 1) {
623 :     fprintf(fh, "after-call: ");
624 :     display_hits(fh);
625 :     }
626 :     /* once we have decided to call a region, we take the kmers for fI and
627 :     add them to the counts maintained to assign an OTU to the sequence */
628 :     for (i=0; (i <= last_hit); i++) {
629 :     if (hits[i].fI == current_fI) {
630 :     int j;
631 :     for (j=0; (j < num_oI) && (oI_counts[j].oI != hits[i].oI); j++) {}
632 :     if (j == num_oI) {
633 :     if (num_oI == OI_BUFSZ) {
634 :     j--; /* we overwrite the last entry */
635 :     }
636 :     else
637 :     num_oI++;
638 :     oI_counts[j].oI = hits[i].oI;
639 :     oI_counts[j].count = 1;
640 :     }
641 :     else {
642 :     oI_counts[j].count++;
643 :     }
644 :     /* now we bubble the count back, allowing it to establish */
645 :     while ((j > 0) && (oI_counts[j-1].count <= oI_counts[j].count)) {
646 :     int oI_tmp = oI_counts[j-1].oI;
647 :     int count_tmp = oI_counts[j-1].count;
648 :     oI_counts[j-1].oI = oI_counts[j].oI;
649 :     oI_counts[j-1].count = oI_counts[j].count;
650 :     oI_counts[j].oI = oI_tmp;
651 :     oI_counts[j].count = count_tmp;
652 :     j--;
653 :     }
654 :     }
655 :     }
656 :     }
657 :    
658 :     if ((hits[num_hits-2].fI != current_fI) && (hits[num_hits-2].fI == hits[num_hits-1].fI)) {
659 :     current_fI = hits[num_hits-1].fI;
660 :     /* now copy the last two entries to the start of the hits array. Sorry this is so clumsy */
661 :     hits[0].oI = hits[num_hits-2].oI;
662 :     hits[0].from0_in_prot = hits[num_hits-2].from0_in_prot;
663 :     hits[0].avg_off_from_end = hits[num_hits-2].avg_off_from_end;
664 :     hits[0].fI = hits[num_hits-2].fI;
665 :     hits[0].function_wt = hits[num_hits-2].function_wt;
666 :    
667 :     hits[1].oI = hits[num_hits-1].oI;
668 :     hits[1].from0_in_prot = hits[num_hits-1].from0_in_prot;
669 :     hits[1].avg_off_from_end = hits[num_hits-1].avg_off_from_end;
670 :     hits[1].fI = hits[num_hits-1].fI;
671 :     hits[1].function_wt = hits[num_hits-1].function_wt;
672 :     num_hits = 2;
673 :     }
674 :     else {
675 :     num_hits = 0;
676 :     }
677 :     }
678 :    
679 :     void KmerGuts::gather_hits(int ln_DNA, char strand,int prot_off,char *pseq,
680 :     unsigned char *pIseq, kmer_handle_t *kmersH, FILE *fh) {
681 :    
682 :     if (debug >= 3) {
683 :     fprintf(fh, "translated: %c\t%d\t%s\n",strand,prot_off,pseq);
684 :     }
685 :    
686 :     unsigned char *p = pIseq;
687 :     /* pseq and pIseq are the same length */
688 :    
689 :     unsigned char *bound = pIseq + strlen(pseq) - K;
690 :     advance_past_ambig(&p,bound);
691 :     unsigned long long encodedK=0;
692 :     if (p < bound) {
693 :     encodedK = encoded_kmer(p);
694 :     }
695 :     while (p < bound) {
696 :     long long where = lookup_hash_entry(kmersH->kmer_table,encodedK);
697 :     if (where >= 0) {
698 :     sig_kmer_t *kmers_hash_entry = &(kmersH->kmer_table[where]);
699 :     int avg_off_end = kmers_hash_entry->avg_from_end;
700 :     int fI = kmers_hash_entry->function_index;
701 :     int oI = kmers_hash_entry->otu_index;
702 :     float f_wt = kmers_hash_entry->function_wt;
703 :     if (debug >= 1) {
704 :     if (hits_only)
705 :     fprintf(fh, "%ld\t%s\n",encodedK, current_id);
706 :     else
707 :     fprintf(fh, "HIT\t%ld\t%lld\t%d\t%d\t%0.3f\t%d\n",p-pIseq,encodedK,avg_off_end,fI,f_wt,oI);
708 :     }
709 :    
710 :     if ((num_hits > 0) && (hits[num_hits-1].from0_in_prot + max_gap) < (p-pIseq)) {
711 :     if (num_hits >= min_hits) {
712 :     process_set_of_hits(kmersH, fh);
713 :     }
714 :     else {
715 :     num_hits = 0;
716 :     }
717 :     }
718 :    
719 :     if (num_hits == 0) {
720 :     current_fI = fI; /* if this is the first, set the current_fI */
721 :     }
722 :    
723 :     if ((! order_constraint) || (num_hits == 0) ||
724 :     ((fI == hits[num_hits-1].fI) &&
725 :     (abs(((p-pIseq) - hits[num_hits-1].from0_in_prot) -
726 :     (hits[num_hits-1].avg_off_from_end - avg_off_end)
727 :     ) <= 20))) {
728 :     /* we have a new hit, so we add it to the global set of hits */
729 :     hits[num_hits].oI = oI;
730 :     hits[num_hits].fI = fI;
731 :     hits[num_hits].from0_in_prot = p-pIseq;
732 :     hits[num_hits].avg_off_from_end = avg_off_end;
733 :     hits[num_hits].function_wt = f_wt;
734 :     if (num_hits < MAX_HITS_PER_SEQ - 2)
735 :     num_hits++;
736 :     if (debug > 1) {
737 :     fprintf(fh, "after-hit: ");
738 :     display_hits(fh);
739 :     }
740 :     if ((num_hits > 1) && (current_fI != fI) && /* if we have a pair of new fIs, it is time to */
741 :     (hits[num_hits-2].fI == hits[num_hits-1].fI)) { /* process one set and initialize the next */
742 :     process_set_of_hits(kmersH, fh);
743 :     }
744 :     }
745 :     }
746 :     p++;
747 :     if (p < bound) {
748 :     if (*(p+K-1) < 20) {
749 :     encodedK = ((encodedK % CORE) * 20L) + *(p+K-1);
750 :     }
751 :     else {
752 :     p += K;
753 :     advance_past_ambig(&p,bound);
754 :     if (p < bound) {
755 :     encodedK = encoded_kmer(p);
756 :     }
757 :     }
758 :     }
759 :     }
760 :     if (num_hits >= min_hits) {
761 :     process_set_of_hits(kmersH, fh);
762 :     }
763 :     num_hits = 0;
764 :     }
765 :    
766 :     void KmerGuts::tabulate_otu_data_for_contig(FILE *fh) {
767 :     int i;
768 :     if (!hits_only)
769 :     {
770 :     fprintf(fh, "OTU-COUNTS\t%s[%d]",current_id,current_length_contig);
771 :     for (i=0; (i < num_oI); i++) {
772 :     fprintf(fh, "\t%d-%d",oI_counts[i].count,oI_counts[i].oI);
773 :     }
774 :     fprintf(fh, "\n");
775 :     }
776 :     num_oI = 0;
777 :     }
778 :    
779 :     void KmerGuts::process_aa_seq(char *id,char *pseq,size_t ln,kmer_handle_t *kmersH, FILE *fh) {
780 :    
781 :     strcpy(current_id,id);
782 :     if (!hits_only)
783 :     fprintf(fh, "PROTEIN-ID\t%s\t%d\n",id,ln);
784 :    
785 :     current_length_contig = ln;
786 :     current_strand = '+';
787 :     current_prot_off = 0;
788 :     int i;
789 :     for (i=0; (i < ln); i++)
790 :     pIseq[i] = to_amino_acid_off(*(pseq+i));
791 :     gather_hits(ln,'+',0,pseq,pIseq,kmersH,fh);
792 :     tabulate_otu_data_for_contig(fh);
793 :     }
794 :    
795 :     void KmerGuts::process_seq(char *id,char *data,kmer_handle_t *kmersH, FILE *fh) {
796 :    
797 :    
798 :     strcpy(current_id,id);
799 :     int ln = strlen(data);
800 :     current_length_contig = ln;
801 :     fprintf(fh, "processing %s[%d]\n",id,ln);
802 :     int i;
803 :     for (i=0; (i < 3); i++) {
804 :    
805 :     translate(data,i,pseq,pIseq);
806 :     current_strand = '+';
807 :     current_prot_off = i;
808 :     if (!hits_only)
809 :     fprintf(fh, "TRANSLATION\t%s\t%d\t%c\t%d\n",current_id,
810 :     current_length_contig,
811 :     current_strand,
812 :     current_prot_off);
813 :     gather_hits(ln,'+',i,pseq,pIseq,kmersH, fh);
814 :     }
815 :     rev_comp(data,cdata);
816 :     for (i=0; (i < 3); i++) {
817 :     translate(cdata,i,pseq,pIseq);
818 :    
819 :     current_strand = '-';
820 :     current_prot_off = i;
821 :     if (!hits_only)
822 :     fprintf(fh, "TRANSLATION\t%s\t%d\t%c\t%d\n",current_id,
823 :     current_length_contig,
824 :     current_strand,
825 :     current_prot_off);
826 :     gather_hits(ln,'-',i,pseq,pIseq,kmersH, fh);
827 :     }
828 :     tabulate_otu_data_for_contig(fh);
829 :     }
830 :    
831 :     void KmerGuts::run_from_filehandle(kmer_handle_t *kmersH, FILE *fh_in, FILE *fh_out)
832 :     {
833 :     int got_gt = 0;
834 :     int i;
835 :     char *p;
836 :     char id[2000];
837 :    
838 :     while (((!got_gt && (fscanf(fh_in,">%s",id) == 1)) ||
839 :     (got_gt && (fscanf(fh_in,"%s",id) == 1)))) {
840 :     while (getc(fh_in) != '\n')
841 :     ;
842 :     if ((id[0] == 'F') &&
843 :     (id[1] == 'L') &&
844 :     (id[2] == 'U') &&
845 :     (id[3] == 'S') &&
846 :     (id[4] == 'H')) { /* ugly compare, since \n is included in id */
847 :    
848 :     /* int rc = fflush(fh_in);
849 :     if (rc != 0) {
850 :     fprintf(stderr,"fflush did not seem to work\n");
851 :     }
852 :     */
853 :     fprintf(fh_out, "//\n");
854 :     got_gt = 0;
855 :     }
856 :     else {
857 :    
858 :     for (p=data; ((i = getc(fh_in)) != -1) && (i != '>');) {
859 :     if ((i != ' ') && (i != '\n'))
860 :     *(p++) = toupper(i);
861 :     }
862 :    
863 :     if (i == '>')
864 :     got_gt = 1;
865 :     else
866 :     got_gt = 0;
867 :    
868 :     *p=0;
869 :     size_t len = p - data;
870 :     if ((len) > MAX_SEQ_LEN) {
871 :     fprintf(stderr,"The contig size exceeds %d; bump MAX_SEQ_LEN\n",MAX_SEQ_LEN);
872 :     exit(1);
873 :     }
874 :    
875 :     /* fprintf(stderr,"%d bytes read\nends with %s\n",(p-data),p-50); */
876 :    
877 :     if (! aa)
878 :     process_seq(id,data,kmersH, fh_out);
879 :     else
880 :     process_aa_seq(id,data, len, kmersH, fh_out);
881 :     fflush(fh_out);
882 :     }
883 :     }
884 :    
885 :     if (debug >= 2)
886 :     fprintf(fh_out, "tot_lookups=%d retry=%d\n",tot_lookups,retry);
887 :     }
888 :    
889 :     void KmerGuts::run_accept_loop(kmer_handle_t *kmersH, in_port_t port, char *port_file, pid_t parent)
890 :     {
891 :     int listenfd = 0, connfd = 0;
892 :     struct sockaddr_in serv_addr;
893 :    
894 :     signal(SIGPIPE, SIG_IGN);
895 :    
896 :     listenfd = socket(AF_INET, SOCK_STREAM, 0);
897 :     memset(&serv_addr, 0, sizeof(serv_addr));
898 :    
899 :     serv_addr.sin_family = AF_INET;
900 :     serv_addr.sin_addr.s_addr = htonl(INADDR_ANY);
901 :     serv_addr.sin_port = htons(port);
902 :    
903 :     if (bind(listenfd, (struct sockaddr*)&serv_addr, sizeof(serv_addr)) < 0)
904 :     {
905 :     perror("bind failed");
906 :     exit(1);
907 :     }
908 :    
909 :     struct sockaddr_in my_addr;
910 :     socklen_t my_len = sizeof(my_addr);
911 :     if (getsockname(listenfd, (struct sockaddr *) &my_addr, &my_len) < 0)
912 :     {
913 :     perror("getsockname failed");
914 :     exit(1);
915 :     }
916 :     in_port_t my_port = ntohs(my_addr.sin_port);
917 :     printf("Listening on %d\n", my_port);
918 :     if (port_file[0])
919 :     {
920 :     FILE *fp = fopen(port_file, "w");
921 :     if (!fp)
922 :     {
923 :     fprintf(stderr, "error opening %s for writing: %s\n", port_file, strerror(errno));
924 :     exit(1);
925 :     }
926 :     fprintf(fp, "%d\n", my_port);
927 :     fclose(fp);
928 :     }
929 :    
930 :     listen(listenfd, 10);
931 :    
932 :     int save_aa = aa;
933 :     int save_hits_only = hits_only;
934 :     int save_debug = debug;
935 :     int save_min_hits = min_hits;
936 :     int save_min_weighted_hits = min_weighted_hits;
937 :     int save_order_constraint = order_constraint;
938 :     int save_max_gap = max_gap;
939 :    
940 :     while(1)
941 :     {
942 :     /*
943 :     * If we have a parent set, and parent doesn't exist, exit.
944 :     */
945 :     if (parent > 0)
946 :     {
947 :     int rc = kill(parent, 0);
948 :     if (rc < 0)
949 :     {
950 :     fprintf(stderr, "Parent process %d does not exist any more, exiting\n", parent);
951 :     exit(0);
952 :     }
953 :     }
954 :    
955 :     aa = save_aa;
956 :     hits_only = save_hits_only;
957 :     debug = save_debug;
958 :     min_hits = save_min_hits;
959 :     min_weighted_hits = save_min_weighted_hits;
960 :     order_constraint = save_order_constraint;
961 :     max_gap = save_max_gap;
962 :    
963 :     connfd = accept(listenfd, (struct sockaddr*)NULL, NULL);
964 :     struct sockaddr_in peer;
965 :     socklen_t peer_len = sizeof(peer);
966 :     memset(&peer, 0, sizeof(peer));
967 :    
968 :     getpeername(connfd, (struct sockaddr *) &peer, &peer_len);
969 :    
970 :     char *who = inet_ntoa(peer.sin_addr);
971 :     // fprintf(stderr, "connection from %s\n", who);
972 :    
973 :     FILE *fh_in = fdopen(connfd, "r");
974 :     FILE *fh_out = fdopen(connfd, "w");
975 :    
976 :     /*
977 :     * If the first line starts with a '-', then it's
978 :     * a set of options to be set for this document.
979 :     */
980 :    
981 :     int c = fgetc(fh_in);
982 :     if (c == '-')
983 :     {
984 :     char linebuf[1024];
985 :     linebuf[0] = c;
986 :    
987 :     if (fgets(linebuf + 1, sizeof(linebuf) - 2, fh_in) == 0)
988 :     {
989 :     fprintf(stderr, "Error reading options line from %s\n", who);
990 :     fclose(fh_in);
991 :     fclose(fh_out);
992 :     continue;
993 :     }
994 :    
995 :     char *s = strchr(linebuf, '\n');
996 :     if (s)
997 :     *s = 0;
998 :    
999 :     /* parse into an argv */
1000 :     const int max_args = 20;
1001 :     char *argv[max_args + 2];
1002 :     int n = 0;
1003 :     // fprintf(stderr, "Parsing option line '%s'\n", linebuf);
1004 :     argv[n++] = "nothing";
1005 :    
1006 :     s = strtok(linebuf, " \t");
1007 :     while (s)
1008 :     {
1009 :     // fprintf(stderr, "argv[%d] = '%s'\n", n, s);
1010 :     argv[n++] = s;
1011 :     if (n >= max_args)
1012 :     {
1013 :     fprintf(stderr, "too many args in connection from %s\n", who);
1014 :     fprintf(fh_out, "ERR too many args\n");
1015 :     fflush(fh_out);
1016 :     fclose(fh_out);
1017 :     fclose(fh_in);
1018 :     s = 0;
1019 :     break;
1020 :     }
1021 :     s = strtok(0, " \t");
1022 :     }
1023 :     if (n >= max_args)
1024 :     {
1025 :     continue;
1026 :     }
1027 :     argv[n] = 0;
1028 :    
1029 :     char *past;
1030 :     int arg_error = 0;
1031 :    
1032 :     optind = 1;
1033 :     while ((c = getopt(n, argv, "ad:m:M:Og:")) != -1)
1034 :     {
1035 :     switch (c) {
1036 :     case 'a':
1037 :     aa = 1;
1038 :     break;
1039 :     case 'd':
1040 :     debug = strtol(optarg,&past,0);
1041 :     break;
1042 :     case 'm':
1043 :     min_hits = strtol(optarg,&past,0);
1044 :     break;
1045 :     case 'M':
1046 :     min_weighted_hits = strtol(optarg,&past,0);
1047 :     break;
1048 :     case 'O':
1049 :     order_constraint = 1;
1050 :     break;
1051 :     case 'g':
1052 :     max_gap = strtol(optarg,&past,0);
1053 :     break;
1054 :     default:
1055 :     fprintf(fh_out, "ERR invalid argument %c\n", c);
1056 :     fflush(fh_out);
1057 :     fclose(fh_out);
1058 :     fclose(fh_in);
1059 :     arg_error = 1;
1060 :     break;
1061 :     }
1062 :     }
1063 :     if (arg_error)
1064 :     continue;
1065 :     if (!hits_only)
1066 :     fprintf(fh_out, "OK aa=%d debug=%d min_hits=%d min_weighted_hits=%d order_constraint=%d max_gap=%d\n",
1067 :     aa, debug, min_hits, min_weighted_hits, order_constraint, max_gap);
1068 :     }
1069 :     else
1070 :     {
1071 :     ungetc(c, fh_in);
1072 :     }
1073 :    
1074 :    
1075 :     run_from_filehandle(kmersH, fh_in, fh_out);
1076 :    
1077 :     fflush(fh_out);
1078 :     fclose(fh_in);
1079 :     fclose(fh_out);
1080 :    
1081 :     }
1082 :     }

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3