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

Annotation of /Kmers2/rna_seq.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.1 - (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 :     char *vals = malloc(ID_TAB_SZ);
327 :     char *p = vals;
328 :     FILE *ifp = fopen(filename,"r");
329 :     if (ifp == NULL) {
330 :     printf("could not open %s\n",filename);
331 :     exit(1);
332 :     }
333 :    
334 :     int sz = 0;
335 :     int j;
336 :     while ((fscanf(ifp,"%d\t",&j) == 1) && fgets(p,1000,ifp)) {
337 :     if (sz != j) {
338 :     fprintf(stderr,"Your index must be dense and in order (see line %ld, should be %d)\n",p-vals,sz);
339 :     exit(1);
340 :     }
341 :     *(id_index_ar[sz]) = p; /* the fgets leaves the \n at the end of each line */
342 :     p += strlen(*(id_index_ar[sz])) -1;
343 :     *(p++) = '\0';
344 :     if ((sz >= ID_TAB_MAX_ENTRIES) || ((p-vals) > (ID_TAB_SZ - 1000))) {
345 :     fprintf(stderr,"Your contig id array is too small; bump ID_TAB_MAX_ENTRIES and ID_TAB_SZ\n");
346 :     exit(1);
347 :     }
348 :     sz += 1;
349 :     }
350 :    
351 :     int c;
352 :     while ((c = getc(ifp)) != '\n') {}
353 :    
354 :     char kmer_string[K+1];
355 :     int id;
356 :     int pos;
357 :     while (fscanf(ifp,"%s\t%d\t%d",kmer_string,&id,&pos) == 3) {
358 :     raw_t v = encode_plus(kmer_string);
359 :     if (! v.bad_input) {
360 :     hash_entry_t *hpos = *hash_tab + v.h1;
361 :     if (hpos->kmer_position == 0) {
362 :     hpos->id_index = id;
363 :     hpos->kmer_position = pos;
364 :     hpos->kmer_hash2 = v.h2;
365 :     }
366 :     else {
367 :     /* fprintf(stderr,"clash: losing %s\n",kmer_string); */
368 :     }
369 :     v = encode_minus(kmer_string);
370 :     hpos = *hash_tab + v.h1;
371 :     if (hpos->kmer_position == 0) {
372 :     hpos->id_index = id;
373 :     hpos->kmer_position = -pos;
374 :     hpos->kmer_hash2 = v.h2;
375 :     }
376 :     else {
377 :     /* fprintf(stderr,"clash: losing comlement of %s\n",kmer_string); */
378 :     }
379 :     }
380 :     }
381 :     }
382 :    
383 :     buffer_t *buffer_loc;
384 :    
385 :     void allocate_buffer() {
386 :    
387 :     size_t sz = sizeof(buffer_t);
388 :     buffer_loc = (buffer_t *) malloc(sz);
389 :     }
390 :    
391 :     void print_usage();
392 :    
393 :     char *kmers;
394 :     char *input;
395 :     char *output;
396 :    
397 :     hash_entry_t *ref_hash;
398 :     char **ref_id_index_ar;
399 :    
400 :     int read_input(handle_t *h) {
401 :     buffer_t *bp = buffer_loc;
402 :     int fastaII = 0;
403 :     bp->data_in_buffer = 0;
404 :     char *next_id = bp->ids_in_buffer;
405 :    
406 :     int got = 0;
407 :     while ((bp->data_in_buffer < (WU_PER_BUFFER - MAX_READ_SIZE)) &&
408 :     get_id(h,next_id) &&
409 :     (get_seq(h,bp->seq+bp->data_in_buffer,
410 :     (WU_PER_BUFFER+(K-1) - bp->data_in_buffer),
411 :     &(bp->where_seq[fastaII].ln)))) {
412 :    
413 :     got++;
414 :     bp->where_seq[fastaII].id = next_id;
415 :     bp->where_seq[fastaII].fastaB = 0;
416 :     bp->where_seq[fastaII].bufferB = bp->data_in_buffer;
417 :     bp->where_seq[fastaII+1].fastaB = -1;
418 :     next_id += strlen(next_id)+1;
419 :     bp->data_in_buffer += bp->where_seq[fastaII].ln;
420 :     if (bp->where_seq[fastaII].ln > K) {
421 :     bp->where_seq[fastaII].ln -= (K-1);
422 :     }
423 :     else {
424 :     fprintf(stderr,"BUG\n");
425 :     }
426 :    
427 :     fastaII++;
428 :     bp->where_seq[fastaII].fastaB = -1;
429 :     }
430 :     return got;
431 :     }
432 :    
433 :     void work() {
434 :     buffer_t *bp = buffer_loc;
435 :     int i;
436 :     for (i=0; (i < bp->data_in_buffer);i++) {
437 :     bp->raw[i] = encode_plus(bp->seq+i);
438 :     }
439 :     }
440 :    
441 :     /* ============== GLOBAL: hit recording =================== */
442 :     /* We consider a string of hits to begin with two against the same
443 :     contig, with a close offset.
444 :     */
445 :    
446 :     char id_of_read[300];
447 :    
448 :     typedef struct hit {
449 :     int contig_id_index;
450 :     int contig_offset;
451 :     int strand;
452 :     char * read_id;
453 :     int read_offset;
454 :     } hit_t;
455 :    
456 :     hit_t first;
457 :     hit_t last;
458 :     hit_t possible_challenger;
459 :    
460 :     int num_hits;
461 :    
462 :     /* ========= end hit global ========== */
463 :    
464 :     void init_hit_list(char *read_id) {
465 :     if (read_id == NULL) { fprintf(stderr,"BAD\n"); }
466 :     strcpy(id_of_read,read_id);
467 :     possible_challenger.contig_id_index = -1;
468 :     num_hits = 0;
469 :     }
470 :    
471 :     void flush_hit_list() {
472 :     if (num_hits > 1) {
473 :     if (first.contig_id_index == -1) {
474 :     fprintf(stderr,"BAD BAD\n"); exit(1);
475 :     }
476 :     if (first.strand == 0) {
477 :     fprintf(stdout,"%d\t%d\t%s\t%d\t%d\t%d\t%d\t%d\n",
478 :     num_hits,
479 :     first.contig_id_index,
480 :     first.read_id,
481 :     first.strand,
482 :     first.contig_offset+1,
483 :     last.contig_offset+K,
484 :     first.read_offset+1,
485 :     last.read_offset+K);
486 :     }
487 :     else {
488 :     fprintf(stdout,"%d\t%d\t%s\t%d\t%d\t%d\t%d\t%d\n",
489 :     num_hits,
490 :     first.contig_id_index,
491 :     first.read_id,
492 :     first.strand,
493 :     first.contig_offset+((K-1)+1),
494 :     last.contig_offset+1,
495 :     first.read_offset+1,
496 :     last.read_offset+K);
497 :     }
498 :     num_hits = 0;
499 :     }
500 :     }
501 :    
502 :     void record_hit(int off_into_read,int strand,int contig_idI,int off_into_contig) {
503 :     /* fprintf(stderr,"HIT: %s %d %d %d\n",id_of_read,off_into_read,strand,off_into_contig);*/
504 :     if (num_hits == 0) {
505 :     first.contig_id_index = contig_idI;
506 :     first.contig_offset = last.contig_offset = off_into_contig;
507 :     first.strand = strand;
508 :     first.read_id = id_of_read;
509 :     first.read_offset = off_into_read;
510 :     num_hits = 1;
511 :     }
512 :     else {
513 :     int read_gap = abs(off_into_read - first.read_offset);
514 :     if ((first.contig_id_index == contig_idI) &&
515 :     (first.strand == strand) &&
516 :     (abs(last.contig_offset - off_into_contig) < (read_gap + 20)) &&
517 :     (abs(first.contig_offset - off_into_contig) < (read_gap + 20)) &&
518 :     (first.read_id == id_of_read)) {
519 :     last.contig_offset = off_into_contig;
520 :     last.read_offset = off_into_read;
521 :     num_hits++;
522 :     possible_challenger.contig_id_index = -1;
523 :     }
524 :     else {
525 :     if (possible_challenger.contig_id_index == -1) {
526 :     possible_challenger.contig_id_index = contig_idI;
527 :     possible_challenger.contig_offset = off_into_contig;
528 :     possible_challenger.strand = strand;
529 :     possible_challenger.read_id = id_of_read;
530 :     possible_challenger.read_offset = off_into_read;
531 :     }
532 :     else {
533 :     flush_hit_list();
534 :     first = possible_challenger;
535 :     num_hits = 1;
536 :     possible_challenger.contig_id_index = -1;
537 :     record_hit(off_into_read,strand,contig_idI,off_into_contig);
538 :     }
539 :     }
540 :     }
541 :     }
542 :    
543 :     void process_raw(int p1, int p2, int fasta_off) {
544 :     hash_entry_t *t = ref_hash+p1;
545 :    
546 :     if (p2 == t->kmer_hash2) {
547 :     if (t->kmer_position < 0) {
548 :     record_hit(fasta_off,1,t->id_index,-(t->kmer_position));
549 :     }
550 :     else {
551 :     record_hit(fasta_off,0,t->id_index,t->kmer_position);
552 :     }
553 :     }
554 :     }
555 :    
556 :     void write_output() {
557 :    
558 :     buffer_t *bp = buffer_loc;
559 :     fasta_mapping_t *fm = bp->where_seq;
560 :     while (fm->fastaB != -1) {
561 :     init_hit_list(fm->id);
562 :     raw_t *outP = bp->raw+fm->bufferB;
563 :     int fasta_off = fm->fastaB;
564 :     int entries_to_process = fm->ln;
565 :     while (entries_to_process) {
566 :     if (! outP->bad_input) {
567 :     process_raw(outP->h1,outP->h2,fasta_off); /* strand: 0 -> +, 1 -> - */
568 :     }
569 :     entries_to_process--;
570 :     outP++;
571 :     fasta_off++;
572 :     }
573 :     flush_hit_list();
574 :     fm++;
575 :     }
576 :     }
577 :    
578 :    
579 :     int main(int argc, char *argv[]) {
580 :     int option = 0;
581 :    
582 :     /*
583 :     s - the name of the file containing the contig index table and the
584 :     signature kmers for the reference genome
585 :    
586 :     K - the Kmer size
587 :     */
588 :     while ((option = getopt(argc, argv,"s:K:")) != -1) {
589 :     switch (option) {
590 :     case 's' : kmers = optarg;
591 :     break;
592 :     case 'K' : K = atoi(optarg);
593 :     break;
594 :     default: print_usage();
595 :     exit(EXIT_FAILURE);
596 :     }
597 :     }
598 :    
599 :     allocate_buffer();
600 :     load_kmers(kmers,&ref_hash,&ref_id_index_ar);
601 :     handle_t *h = get_handle();
602 :     while (read_input(h)) {
603 :     work();
604 :     write_output();
605 :     }
606 :     exit(0);
607 :     }
608 :    
609 :     void print_usage() {
610 :     printf("Usage: rna_seq -s signature-kmers -K kmer-size\n");
611 :     }

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3