[Bio] / Kmers2 / rna_seq.c Repository:
ViewVC logotype

Annotation of /Kmers2/rna_seq.c

Parent Directory Parent Directory | Revision Log Revision Log


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

1 : overbeek 1.1 #include <unistd.h>
2 :     #include <stdio.h>
3 :     #include <stdlib.h>
4 :     #include <getopt.h>
5 :     #include <sys/types.h>
6 :     #include <sys/ipc.h>
7 :     #include <stdio.h>
8 :     #include <string.h>
9 :     #include <fcntl.h>
10 :     #include <sys/stat.h>
11 :     #include <sys/wait.h>
12 :     #include <sys/mman.h>
13 :     #include <errno.h> /* errno */
14 :    
15 :     #define READSZ 4*BUFSIZ
16 :    
17 :     int K = 24; /* cannot exceed 39 */
18 :     int NUM_BUFFERS = 20;
19 :    
20 :     #define WU_PER_BUFFER 20000
21 :     #define MAX_READ_SIZE 2000
22 :     #define SEQ_CHUNKS_PER_BUFFER 5000
23 :     #define MAX_SZ_IDS_IN_BUFFER 20000
24 :    
25 :     char trans[256] = {4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,0,4,1,4,4,4,2,4,4,4,4,4,4,4,4,4,4,4,4,3,3,4,4,4,4,4,4,4,4,4,4,4,0,4,1,4,4,4,2,4,4,4,4,4,4,4,4,4,4,4,4,3,3,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4};
26 :    
27 :     char transC[256] = {4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,3,4,2,4,4,4,1,4,4,4,4,4,4,4,4,4,4,4,4,0,0,4,4,4,4,4,4,4,4,4,4,4,3,4,2,4,4,4,1,4,4,4,4,4,4,4,4,4,4,4,4,0,0,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4};
28 :    
29 :     typedef struct raw_output {
30 :     int h1;
31 :     int h2;
32 :     int bad_input;
33 :     } raw_t;
34 :    
35 :     typedef struct fasta_mapping {
36 :     char *id; /* ID array */
37 :     int fastaB; /* -1 -> empty */
38 :     int bufferB;
39 :     int ln;
40 :     } fasta_mapping_t;
41 :    
42 :     typedef struct buffer {
43 :     char seq[WU_PER_BUFFER+40]; /* we include a K-1 character pad */
44 :     int data_in_buffer;
45 :     raw_t raw[WU_PER_BUFFER];
46 :     fasta_mapping_t where_seq[SEQ_CHUNKS_PER_BUFFER];
47 :     char ids_in_buffer[MAX_SZ_IDS_IN_BUFFER];
48 :     } buffer_t;
49 :    
50 :     #define HASH2 15485867
51 :     #define HSIZE 100000007
52 :    
53 :     #define ID_TAB_MAX_ENTRIES 100000
54 :     #define ID_TAB_SZ 1000000
55 :    
56 :     typedef struct hash_entry {
57 :     int kmer_hash2;
58 :     int id_index;
59 :     int kmer_position; /* strand: positive for plus; -(p+(K-1)) for minus */
60 :     } hash_entry_t;
61 :    
62 :     typedef struct handle {
63 :     char buffer[READSZ];
64 :     char *buffP;
65 :     int state;
66 :     int left;
67 :     } handle_t;
68 :    
69 :     int fill (handle_t *h) {
70 :     h->left = read (STDIN_FILENO, h->buffer, READSZ);
71 :     h->buffP = h->buffer;
72 :     return h->left;
73 :     }
74 :    
75 :     handle_t *get_handle() {
76 :     handle_t *h = (handle_t *) malloc(sizeof(handle_t));
77 :    
78 :     char c;
79 :     if ((read(STDIN_FILENO,&c,1) < 1) || (c != '>')) {
80 :     fprintf (stderr, "Input does not start with >\n");
81 :     exit(0);
82 :     }
83 :    
84 :     h->left = 0;
85 :     h->state = 1;
86 :    
87 :     return h;
88 :     }
89 :    
90 :     int get_id(handle_t *h, char *id) {
91 :     char c;
92 :     int got = 0;
93 :    
94 :     while (!got) {
95 :     while (h->left && ((c = *(h->buffP)) != ' ') && (c != '\t') && (c != '\n') ) {
96 :     *id++ = c;
97 :     h->left--;
98 :     h->buffP++;
99 :     }
100 :    
101 :     if (h->left == 0) {
102 :     if (!fill(h)) {
103 :     return 0;
104 :     }
105 :     } else {
106 :     int got2 = 0;
107 :     while (! got2) {
108 :     while (h->left && ((c = *(h->buffP)) != '\n') ) {
109 :     h->left--;
110 :     h->buffP++;
111 :     }
112 :     if (h->left == 0) {
113 :     if (!fill(h)) {
114 :     return 0;
115 :     }
116 :     } else {
117 :     h->left--;
118 :     h->buffP++;
119 :     got2 = 1;
120 :     }
121 :     }
122 :     *id = '\0';
123 :     got = 1;
124 :     }
125 :     }
126 :     return 1;
127 :     }
128 :    
129 :     int get_seq(handle_t *h, char *where, int want, int *how_many) {
130 :    
131 :     *how_many = 0;
132 :     int got = 0;
133 :     char c;
134 :     while (! got) {
135 :     while (want && h->left && ((c = *(h->buffP)) != '>')) {
136 :     if ((c != ' ') && (c != '\t') && (c != '\n')) {
137 :     *where++ = *(h->buffP);
138 :     want--;
139 :     (*how_many)++;
140 :     }
141 :     h->left--;
142 :     h->buffP++;
143 :     }
144 :     if (want == 0) {
145 :     return 0;
146 :     }
147 :     else {
148 :     if (h->left == 0) {
149 :     if (! fill(h)) {
150 :     return 1;
151 :     }
152 :     }
153 :     else {
154 :     h->left--;
155 :     h->buffP++;
156 :     return 1;
157 :     }
158 :     }
159 :     }
160 :     return 0;
161 :     }
162 :    
163 :     int dna_char(c)
164 :     char c;
165 :     {
166 :     switch (c)
167 :     {
168 :     case 'a':
169 :     case 'A':
170 :     return 0;
171 :    
172 :     case 'c':
173 :     case 'C':
174 :     return 1;
175 :    
176 :     case 'g':
177 :     case 'G':
178 :     return 2;
179 :    
180 :     case 't':
181 :     case 'u':
182 :     case 'T':
183 :     case 'U':
184 :     return 3;
185 :    
186 :     default:
187 :     return 4;;
188 :     }
189 :     }
190 :    
191 :     char compl(c)
192 :     char c;
193 :     {
194 :     switch (c)
195 :     {
196 :     case 'a':
197 :     return 't';
198 :     case 'A':
199 :     return 'T';
200 :    
201 :     case 'c':
202 :     return 'g';
203 :     case 'C':
204 :     return 'G';
205 :    
206 :     case 'g':
207 :     return 'c';
208 :     case 'G':
209 :     return 'C';
210 :    
211 :     case 't':
212 :     case 'u':
213 :     return 'a';
214 :     case 'T':
215 :     case 'U':
216 :     return 'A';
217 :    
218 :     case 'm':
219 :     return 'k';
220 :     case 'M':
221 :     return 'K';
222 :    
223 :     case 'r':
224 :     return 'y';
225 :     case 'R':
226 :     return 'Y';
227 :    
228 :     case 'w':
229 :     return 'w';
230 :     case 'W':
231 :     return 'W';
232 :    
233 :     case 's':
234 :     return 'S';
235 :     case 'S':
236 :     return 'S';
237 :    
238 :     case 'y':
239 :     return 'r';
240 :     case 'Y':
241 :     return 'R';
242 :    
243 :     case 'k':
244 :     return 'm';
245 :     case 'K':
246 :     return 'M';
247 :    
248 :     case 'b':
249 :     return 'v';
250 :     case 'B':
251 :     return 'V';
252 :    
253 :     case 'd':
254 :     return 'h';
255 :     case 'D':
256 :     return 'H';
257 :    
258 :     case 'h':
259 :     return 'd';
260 :     case 'H':
261 :     return 'D';
262 :    
263 :     case 'v':
264 :     return 'b';
265 :     case 'V':
266 :     return 'B';
267 :    
268 :     case 'n':
269 :     return 'n';
270 :     case 'N':
271 :     return 'N';
272 :    
273 :     default:
274 :     return c;
275 :     }
276 :     }
277 :    
278 :     raw_t encode_plus(char *kmer) {
279 :     long long v = 0;
280 :     int i;
281 :     raw_t r;
282 :    
283 :     r.bad_input = 0;
284 :     for (i = 0; i < K; i++) {
285 :     int c = trans[*(kmer+i)];
286 :     if (c == 4) {
287 :     r.bad_input = 1;;
288 :     }
289 :     v = (v << 2) | c;
290 :     }
291 :     r.h1 = v % HSIZE;
292 :     r.h2 = v % HASH2;
293 :     return r;
294 :     }
295 :    
296 :     raw_t encode_minus(char *kmer) {
297 :     long long v = 0;
298 :     int i;
299 :     raw_t r;
300 :    
301 :     r.bad_input = 0;
302 :     for (i = K-1; i >= 0; i--) {
303 :     int c = transC[*(kmer+i)];
304 :     if (c == 4) {
305 :     r.bad_input = 1;;
306 :     }
307 :     v = (v << 2) | c;
308 :    
309 :     }
310 :     r.h1 = v % HSIZE;
311 :     r.h2 = v % HASH2;
312 :     return r;
313 :     }
314 :    
315 :     void load_kmers(char *filename, hash_entry_t **hash_tab, char ***id_index_ar) {
316 :    
317 :     *hash_tab = malloc(HSIZE * sizeof(hash_entry_t));
318 :     int i;
319 :     hash_entry_t *p1;
320 :     for (i=0; (i < HSIZE); i++) {
321 :     p1 = *hash_tab + i;
322 :     p1->kmer_position = 0;
323 :     }
324 :    
325 :     *id_index_ar = malloc(ID_TAB_MAX_ENTRIES * sizeof(char *));
326 : overbeek 1.2 char **p_index = *id_index_ar;
327 :    
328 : overbeek 1.1 char *vals = malloc(ID_TAB_SZ);
329 :     char *p = vals;
330 :     FILE *ifp = fopen(filename,"r");
331 :     if (ifp == NULL) {
332 :     printf("could not open %s\n",filename);
333 :     exit(1);
334 :     }
335 :    
336 :     int sz = 0;
337 :     int j;
338 :     while ((fscanf(ifp,"%d\t",&j) == 1) && fgets(p,1000,ifp)) {
339 :     if (sz != j) {
340 :     fprintf(stderr,"Your index must be dense and in order (see line %ld, should be %d)\n",p-vals,sz);
341 :     exit(1);
342 :     }
343 : overbeek 1.2 p_index[sz] = p;
344 :    
345 :     p += strlen(p) -1;
346 : overbeek 1.1 *(p++) = '\0';
347 :     if ((sz >= ID_TAB_MAX_ENTRIES) || ((p-vals) > (ID_TAB_SZ - 1000))) {
348 :     fprintf(stderr,"Your contig id array is too small; bump ID_TAB_MAX_ENTRIES and ID_TAB_SZ\n");
349 :     exit(1);
350 :     }
351 :     sz += 1;
352 :     }
353 :    
354 :     int c;
355 :     while ((c = getc(ifp)) != '\n') {}
356 :    
357 :     char kmer_string[K+1];
358 :     int id;
359 :     int pos;
360 :     while (fscanf(ifp,"%s\t%d\t%d",kmer_string,&id,&pos) == 3) {
361 :     raw_t v = encode_plus(kmer_string);
362 :     if (! v.bad_input) {
363 :     hash_entry_t *hpos = *hash_tab + v.h1;
364 :     if (hpos->kmer_position == 0) {
365 :     hpos->id_index = id;
366 :     hpos->kmer_position = pos;
367 :     hpos->kmer_hash2 = v.h2;
368 :     }
369 :     else {
370 :     /* fprintf(stderr,"clash: losing %s\n",kmer_string); */
371 :     }
372 :     v = encode_minus(kmer_string);
373 :     hpos = *hash_tab + v.h1;
374 :     if (hpos->kmer_position == 0) {
375 :     hpos->id_index = id;
376 :     hpos->kmer_position = -pos;
377 :     hpos->kmer_hash2 = v.h2;
378 :     }
379 :     else {
380 :     /* fprintf(stderr,"clash: losing comlement of %s\n",kmer_string); */
381 :     }
382 :     }
383 :     }
384 :     }
385 :    
386 :     buffer_t *buffer_loc;
387 :    
388 :     void allocate_buffer() {
389 :    
390 :     size_t sz = sizeof(buffer_t);
391 :     buffer_loc = (buffer_t *) malloc(sz);
392 :     }
393 :    
394 :     void print_usage();
395 :    
396 :     char *kmers;
397 :     char *input;
398 :     char *output;
399 :    
400 :     hash_entry_t *ref_hash;
401 :     char **ref_id_index_ar;
402 :    
403 :     int read_input(handle_t *h) {
404 :     buffer_t *bp = buffer_loc;
405 :     int fastaII = 0;
406 :     bp->data_in_buffer = 0;
407 :     char *next_id = bp->ids_in_buffer;
408 :    
409 :     int got = 0;
410 :     while ((bp->data_in_buffer < (WU_PER_BUFFER - MAX_READ_SIZE)) &&
411 :     get_id(h,next_id) &&
412 :     (get_seq(h,bp->seq+bp->data_in_buffer,
413 :     (WU_PER_BUFFER+(K-1) - bp->data_in_buffer),
414 :     &(bp->where_seq[fastaII].ln)))) {
415 :    
416 :     got++;
417 :     bp->where_seq[fastaII].id = next_id;
418 :     bp->where_seq[fastaII].fastaB = 0;
419 :     bp->where_seq[fastaII].bufferB = bp->data_in_buffer;
420 :     bp->where_seq[fastaII+1].fastaB = -1;
421 :     next_id += strlen(next_id)+1;
422 :     bp->data_in_buffer += bp->where_seq[fastaII].ln;
423 :     if (bp->where_seq[fastaII].ln > K) {
424 :     bp->where_seq[fastaII].ln -= (K-1);
425 :     }
426 :     else {
427 :     fprintf(stderr,"BUG\n");
428 :     }
429 :    
430 :     fastaII++;
431 :     bp->where_seq[fastaII].fastaB = -1;
432 :     }
433 :     return got;
434 :     }
435 :    
436 :     void work() {
437 :     buffer_t *bp = buffer_loc;
438 :     int i;
439 :     for (i=0; (i < bp->data_in_buffer);i++) {
440 :     bp->raw[i] = encode_plus(bp->seq+i);
441 :     }
442 :     }
443 :    
444 :     /* ============== GLOBAL: hit recording =================== */
445 :     /* We consider a string of hits to begin with two against the same
446 :     contig, with a close offset.
447 :     */
448 :    
449 :     char id_of_read[300];
450 :    
451 :     typedef struct hit {
452 :     int contig_id_index;
453 :     int contig_offset;
454 :     int strand;
455 :     char * read_id;
456 :     int read_offset;
457 :     } hit_t;
458 :    
459 :     hit_t first;
460 :     hit_t last;
461 :     hit_t possible_challenger;
462 :    
463 :     int num_hits;
464 :    
465 :     /* ========= end hit global ========== */
466 :    
467 :     void init_hit_list(char *read_id) {
468 :     if (read_id == NULL) { fprintf(stderr,"BAD\n"); }
469 :     strcpy(id_of_read,read_id);
470 :     possible_challenger.contig_id_index = -1;
471 :     num_hits = 0;
472 :     }
473 :    
474 :     void flush_hit_list() {
475 :     if (num_hits > 1) {
476 :     if (first.contig_id_index == -1) {
477 :     fprintf(stderr,"BAD BAD\n"); exit(1);
478 :     }
479 :     if (first.strand == 0) {
480 :     fprintf(stdout,"%d\t%d\t%s\t%d\t%d\t%d\t%d\t%d\n",
481 :     num_hits,
482 :     first.contig_id_index,
483 :     first.read_id,
484 :     first.strand,
485 :     first.contig_offset+1,
486 :     last.contig_offset+K,
487 :     first.read_offset+1,
488 :     last.read_offset+K);
489 :     }
490 :     else {
491 :     fprintf(stdout,"%d\t%d\t%s\t%d\t%d\t%d\t%d\t%d\n",
492 :     num_hits,
493 :     first.contig_id_index,
494 :     first.read_id,
495 :     first.strand,
496 :     first.contig_offset+((K-1)+1),
497 :     last.contig_offset+1,
498 :     first.read_offset+1,
499 :     last.read_offset+K);
500 :     }
501 :     num_hits = 0;
502 :     }
503 :     }
504 :    
505 :     void record_hit(int off_into_read,int strand,int contig_idI,int off_into_contig) {
506 :     /* fprintf(stderr,"HIT: %s %d %d %d\n",id_of_read,off_into_read,strand,off_into_contig);*/
507 :     if (num_hits == 0) {
508 :     first.contig_id_index = contig_idI;
509 :     first.contig_offset = last.contig_offset = off_into_contig;
510 :     first.strand = strand;
511 :     first.read_id = id_of_read;
512 :     first.read_offset = off_into_read;
513 :     num_hits = 1;
514 :     }
515 :     else {
516 :     int read_gap = abs(off_into_read - first.read_offset);
517 :     if ((first.contig_id_index == contig_idI) &&
518 :     (first.strand == strand) &&
519 :     (abs(last.contig_offset - off_into_contig) < (read_gap + 20)) &&
520 :     (abs(first.contig_offset - off_into_contig) < (read_gap + 20)) &&
521 :     (first.read_id == id_of_read)) {
522 :     last.contig_offset = off_into_contig;
523 :     last.read_offset = off_into_read;
524 :     num_hits++;
525 :     possible_challenger.contig_id_index = -1;
526 :     }
527 :     else {
528 :     if (possible_challenger.contig_id_index == -1) {
529 :     possible_challenger.contig_id_index = contig_idI;
530 :     possible_challenger.contig_offset = off_into_contig;
531 :     possible_challenger.strand = strand;
532 :     possible_challenger.read_id = id_of_read;
533 :     possible_challenger.read_offset = off_into_read;
534 :     }
535 :     else {
536 :     flush_hit_list();
537 :     first = possible_challenger;
538 :     num_hits = 1;
539 :     possible_challenger.contig_id_index = -1;
540 :     record_hit(off_into_read,strand,contig_idI,off_into_contig);
541 :     }
542 :     }
543 :     }
544 :     }
545 :    
546 :     void process_raw(int p1, int p2, int fasta_off) {
547 :     hash_entry_t *t = ref_hash+p1;
548 :    
549 :     if (p2 == t->kmer_hash2) {
550 :     if (t->kmer_position < 0) {
551 :     record_hit(fasta_off,1,t->id_index,-(t->kmer_position));
552 :     }
553 :     else {
554 :     record_hit(fasta_off,0,t->id_index,t->kmer_position);
555 :     }
556 :     }
557 :     }
558 :    
559 :     void write_output() {
560 :    
561 :     buffer_t *bp = buffer_loc;
562 :     fasta_mapping_t *fm = bp->where_seq;
563 :     while (fm->fastaB != -1) {
564 :     init_hit_list(fm->id);
565 :     raw_t *outP = bp->raw+fm->bufferB;
566 :     int fasta_off = fm->fastaB;
567 :     int entries_to_process = fm->ln;
568 :     while (entries_to_process) {
569 :     if (! outP->bad_input) {
570 :     process_raw(outP->h1,outP->h2,fasta_off); /* strand: 0 -> +, 1 -> - */
571 :     }
572 :     entries_to_process--;
573 :     outP++;
574 :     fasta_off++;
575 :     }
576 :     flush_hit_list();
577 :     fm++;
578 :     }
579 :     }
580 :    
581 :    
582 :     int main(int argc, char *argv[]) {
583 :     int option = 0;
584 :    
585 :     /*
586 :     s - the name of the file containing the contig index table and the
587 :     signature kmers for the reference genome
588 :    
589 :     K - the Kmer size
590 :     */
591 :     while ((option = getopt(argc, argv,"s:K:")) != -1) {
592 :     switch (option) {
593 :     case 's' : kmers = optarg;
594 :     break;
595 :     case 'K' : K = atoi(optarg);
596 :     break;
597 :     default: print_usage();
598 :     exit(EXIT_FAILURE);
599 :     }
600 :     }
601 :    
602 :     allocate_buffer();
603 :     load_kmers(kmers,&ref_hash,&ref_id_index_ar);
604 :     handle_t *h = get_handle();
605 :     while (read_input(h)) {
606 :     work();
607 :     write_output();
608 :     }
609 :     exit(0);
610 :     }
611 :    
612 :     void print_usage() {
613 :     printf("Usage: rna_seq -s signature-kmers -K kmer-size\n");
614 :     }

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3