1 /* read_dna.c 2012-12-27 */
3 /* Copyright 2008-2012 Emmanuel Paradis */
5 /* This file is part of the R-package `ape'. */
6 /* See the file ../COPYING for licensing issues. */
9 #include <Rinternals.h>
11 // The initial code defining and initialising the translation table:
13 // for (i = 0; i < 122; i++) tab_trans[i] = 0x00;
15 // tab_trans[65] = 0x88; /* A */
16 // tab_trans[71] = 0x48; /* G */
17 // tab_trans[67] = 0x28; /* C */
18 // tab_trans[84] = 0x18; /* T */
19 // tab_trans[82] = 0xc0; /* R */
20 // tab_trans[77] = 0xa0; /* M */
21 // tab_trans[87] = 0x90; /* W */
22 // tab_trans[83] = 0x60; /* S */
23 // tab_trans[75] = 0x50; /* K */
24 // tab_trans[89] = 0x30; /* Y */
25 // tab_trans[86] = 0xe0; /* V */
26 // tab_trans[72] = 0xb0; /* H */
27 // tab_trans[68] = 0xd0; /* D */
28 // tab_trans[66] = 0x70; /* B */
29 // tab_trans[78] = 0xf0; /* N */
31 // tab_trans[97] = 0x88; /* a */
32 // tab_trans[103] = 0x48; /* g */
33 // tab_trans[99] = 0x28; /* c */
34 // tab_trans[116] = 0x18; /* t */
35 // tab_trans[114] = 0xc0; /* r */
36 // tab_trans[109] = 0xa0; /* m */
37 // tab_trans[119] = 0x90; /* w */
38 // tab_trans[115] = 0x60; /* s */
39 // tab_trans[107] = 0x50; /* k */
40 // tab_trans[121] = 0x30; /* y */
41 // tab_trans[118] = 0xe0; /* v */
42 // tab_trans[104] = 0xb0; /* h */
43 // tab_trans[100] = 0xd0; /* d */
44 // tab_trans[98] = 0x70; /* b */
45 // tab_trans[110] = 0xf0; /* n */
47 // tab_trans[45] = 0x04; /* - */
48 // tab_trans[63] = 0x02; /* ? */
50 static const unsigned char tab_trans[] = {
51 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, /* 0-9 */
52 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, /* 10-19 */
53 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, /* 20-29 */
54 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, /* 30-39 */
55 0x00, 0x00, 0x00, 0x00, 0x00, 0x04, 0x00, 0x00, 0x00, 0x00, /* 40-49 */
56 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, /* 50-59 */
57 0x00, 0x00, 0x00, 0x02, 0x00, 0x88, 0x70, 0x28, 0xd0, 0x00, /* 60-69 */
58 0x48, 0x00, 0xb0, 0x00, 0x00, 0x50, 0x00, 0xa0, 0xf0, 0x00, /* 70-79 */
59 0x00, 0x00, 0xc0, 0x60, 0x18, 0x00, 0xe0, 0x90, 0x00, 0x30, /* 80-89 */
60 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x88, 0x70, 0x28, /* 90-99 */
61 0xd0, 0x00, 0x00, 0x48, 0xb0, 0x00, 0x00, 0x50, 0x00, 0xa0, /* 100-109 */
62 0xf0, 0x00, 0x00, 0x00, 0xc0, 0x60, 0x18, 0x00, 0xe0, 0x90, /* 110-119 */
63 0x00, 0x30, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, /* 120-129 */
64 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, /* 130-139 */
65 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, /* 140-149 */
66 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, /* 150-159 */
67 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, /* 160-169 */
68 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, /* 170-179 */
69 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, /* 180-189 */
70 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, /* 190-199 */
71 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, /* 200-209 */
72 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, /* 210-219 */
73 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, /* 220-229 */
74 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, /* 230-239 */
75 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, /* 240-249 */
76 0x00, 0x00, 0x00, 0x00, 0x00, 0x00}; /* 250-255 */
78 static const unsigned char hook = 0x3e;
79 static const unsigned char lineFeed = 0x0a;
80 /* static const unsigned char space = 0x20; */
82 SEXP rawStreamToDNAbin(SEXP x)
85 unsigned char *xr, *rseq, *buffer, tmp;
88 PROTECT(x = coerceVector(x, RAWSXP));
92 /* do a 1st pass to find the number of sequences */
95 j = 0; /* use j as a flag */
96 for (i = 0; i < N; i++) {
97 if (j && xr[i] == lineFeed) {
100 } else if (xr[i] == hook) j = 1;
103 PROTECT(obj = allocVector(VECSXP, n));
104 PROTECT(nms = allocVector(STRSXP, n));
106 /* Refine the way the size of the buffer is set? */
107 buffer = (unsigned char *)R_alloc(N, sizeof(unsigned char *));
111 /* 1st read the label... */
115 while (xr[i] != lineFeed) buffer[k++] = xr[i++];
117 SET_STRING_ELT(nms, j, mkChar(buffer));
119 /* ... then read the sequence */
121 while (xr[i] != hook && i < N) {
122 tmp = tab_trans[xr[i++]];
123 /* If we are sure that the FASTA file is correct (ie, the sequence on
124 a single line and only the IUAPC code (plus '-' and '?') is used,
125 then the following check would not be needed; additionally the size
126 of tab_trans could be restriced to 0-121. This check has the
127 advantage that all invalid characters are simply ignored without
129 if (tmp) buffer[n++] = tmp;
131 PROTECT(seq = allocVector(RAWSXP, n));
133 for (k = 0; k < n; k++) rseq[k] = buffer[k];
134 SET_VECTOR_ELT(obj, j, seq);
138 setAttrib(obj, R_NamesSymbol, nms);