src/read_dna.c

   1 /* read_dna.c       2013-01-04 */
   2
   3 /* Copyright 2013 Emmanuel Paradis */
   4
   5 /* This file is part of the R-package `ape'. */
   6 /* See the file ../COPYING for licensing issues. */
   7
   8 #include <R.h>
   9 #include <Rinternals.h>
  10
  11 // The initial code defining and initialising the translation table:
  12 //
  13 //      for (i = 0; i < 122; i++) tab_trans[i] = 0x00;
  14 //
  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 */
  30 //
  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 */
  46 //
  47 //      tab_trans[45] = 0x04; /* - */
  48 //      tab_trans[63] = 0x02; /* ? */
  49
  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 */
  77
  78 static const unsigned char hook = 0x3e;
  79 static const unsigned char lineFeed = 0x0a;
  80 /* static const unsigned char space = 0x20; */
  81
  82 SEXP rawStreamToDNAbin(SEXP x)
  83 {
  84         int N, i, j, k, n, startOfSeq;
  85         unsigned char *xr, *rseq, *buffer, tmp;
  86         SEXP obj, nms, seq;
  87
  88         PROTECT(x = coerceVector(x, RAWSXP));
  89         N = LENGTH(x);
  90         xr = RAW(x);
  91
  92 /* do a 1st pass to find the number of sequences
  93
  94    this code should be robust to '>' present inside
  95    a label or in the header text before the sequences */
  96
  97         n = j = 0; /* use j as a flag */
  98         if (xr[0] == hook) {
  99                 j = 1;
 100                 startOfSeq = 0;
 101         }
 102         i = 1;
 103         for (i = 1; i < N; i++) {
 104                 if (j && xr[i] == lineFeed) {
 105                         n++;
 106                         j = 0;
 107                 } else if (xr[i] == hook) {
 108                         if (!n) startOfSeq = i;
 109                         j = 1;
 110                 }
 111         }
 112
 113         PROTECT(obj = allocVector(VECSXP, n));
 114         PROTECT(nms = allocVector(STRSXP, n));
 115
 116 /* Refine the way the size of the buffer is set? */
 117         buffer = (unsigned char *)R_alloc(N, sizeof(unsigned char *));
 118
 119         i = startOfSeq;
 120         j = 0; /* gives the index of the sequence */
 121         while (i < N) {
 122                 /* 1st read the label... */
 123                 i++;
 124                 k = 0;
 125                 while (xr[i] != lineFeed) buffer[k++] = xr[i++];
 126                 buffer[k] = '\0';
 127                 SET_STRING_ELT(nms, j, mkChar(buffer));
 128                 /* ... then read the sequence */
 129                 n = 0;
 130                 while (xr[i] != hook && i < N) {
 131                         tmp = tab_trans[xr[i++]];
 132 /* If we are sure that the FASTA file is correct (ie, the sequence on
 133    a single line and only the IUAPC code (plus '-' and '?') is used,
 134    then the following check would not be needed; additionally the size
 135    of tab_trans could be restriced to 0-121. This check has the
 136    advantage that all invalid characters are simply ignored without
 137    causing error -- except if '>' occurs in the middle of a sequence. */
 138                         if (tmp) buffer[n++] = tmp;
 139                 }
 140                 PROTECT(seq = allocVector(RAWSXP, n));
 141                 rseq = RAW(seq);
 142                 for (k = 0; k < n; k++) rseq[k] = buffer[k];
 143                 SET_VECTOR_ELT(obj, j, seq);
 144                 UNPROTECT(1);
 145                 j++;
 146         }
 147         setAttrib(obj, R_NamesSymbol, nms);
 148         UNPROTECT(3);
 149         return obj;
 150 }