src/parse.c

   1
   2 #include <stdbool.h>
   3 #include <stdio.h>
   4 #include <string.h>
   5
   6 #include "parse.h"
   7 #include "common.h"
   8
   9
  10 static void str_init(str_t* str)
  11 {
  12     str->n = 0;
  13     str->size = 128;
  14     str->s = malloc_or_die(str->size = 0);
  15     str->s[0] = '\0';
  16 }
  17
  18
  19 static void str_free(str_t* str)
  20 {
  21     free(str->s);
  22 }
  23
  24
  25 /* Reserve space for `size` more characters. */
  26 static void str_reserve_extra(str_t* str, size_t size)
  27 {
  28     if (str->n + size > str->size) {
  29         if (str->n + size > 2 * str->size) {
  30             str->size = str->n + size;
  31         }
  32         else str->size *= 2;
  33         str->s = realloc_or_die(str->s, str->size * sizeof(char));
  34     }
  35 }
  36
  37
  38 /* Copy n characters from c to the end of str. */
  39 static void str_append(str_t* str, char* c, size_t n)
  40 {
  41     str_reserve_extra(str, n);
  42     memcpy(str->s + str->n, c, n);
  43     str->n += n;
  44     str->s[str->n] = '\0';
  45 }
  46
  47
  48 seq_t* seq_create()
  49 {
  50     seq_t* seq = malloc_or_die(sizeof(seq_t));
  51     str_init(&seq->id1);
  52     str_init(&seq->seq);
  53     str_init(&seq->id2);
  54     str_init(&seq->qual);
  55     return seq;
  56 }
  57
  58
  59 void seq_free(seq_t* seq)
  60 {
  61     str_free(&seq->id1);
  62     str_free(&seq->seq);
  63     str_free(&seq->id2);
  64     str_free(&seq->qual);
  65     free(seq);
  66 }
  67
  68
  69 /* This is MurmurHash3. The original C++ code was placed in the public domain
  70  * by its author, Austin Appleby. */
  71
  72 static inline uint32_t fmix(uint32_t h)
  73 {
  74     h ^= h >> 16;
  75     h *= 0x85ebca6b;
  76     h ^= h >> 13;
  77     h *= 0xc2b2ae35;
  78     h ^= h >> 16;
  79
  80     return h;
  81 }
  82
  83
  84 static inline uint32_t rotl32(uint32_t x, int8_t r)
  85 {
  86     return (x << r) | (x >> (32 - r));
  87 }
  88
  89
  90 uint32_t murmurhash3(const uint8_t* data, size_t len_)
  91 {
  92     const int len = (int) len_;
  93     const int nblocks = len / 4;
  94
  95     uint32_t h1 = 0xc062fb4a;
  96
  97     uint32_t c1 = 0xcc9e2d51;
  98     uint32_t c2 = 0x1b873593;
  99
 100     //----------
 101     // body
 102
 103     const uint32_t * blocks = (const uint32_t*) (data + nblocks * 4);
 104
 105     int i;
 106     for(i = -nblocks; i; i++)
 107     {
 108         uint32_t k1 = blocks[i];
 109
 110         k1 *= c1;
 111         k1 = rotl32(k1, 15);
 112         k1 *= c2;
 113
 114         h1 ^= k1;
 115         h1 = rotl32(h1, 13);
 116         h1 = h1*5+0xe6546b64;
 117     }
 118
 119     //----------
 120     // tail
 121
 122     const uint8_t * tail = (const uint8_t*)(data + nblocks*4);
 123
 124     uint32_t k1 = 0;
 125
 126     switch(len & 3)
 127     {
 128         case 3: k1 ^= tail[2] << 16;
 129         case 2: k1 ^= tail[1] << 8;
 130         case 1: k1 ^= tail[0];
 131               k1 *= c1; k1 = rotl32(k1,15); k1 *= c2; h1 ^= k1;
 132     }
 133
 134     //----------
 135     // finalization
 136
 137     h1 ^= len;
 138
 139     h1 = fmix(h1);
 140
 141     return h1;
 142 }
 143
 144
 145 uint32_t seq_hash(const seq_t* seq)
 146 {
 147     /* TODO */
 148     return 0;
 149 }
 150
 151 static const size_t parser_buf_size = 1000000;
 152
 153
 154 struct fastq_t_
 155 {
 156     FILE* file;
 157     size_t readlen;
 158     char* buf;
 159     char* next;
 160     bool linestart;
 161 };
 162
 163
 164
 165 fastq_t* fastq_create(FILE* file)
 166 {
 167     fastq_t* f = malloc_or_die(sizeof(fastq_t));
 168     f->file = file;
 169     f->next = f->buf = malloc_or_die(parser_buf_size);
 170     f->readlen = 0;
 171     f->linestart = true;
 172     return f;
 173 }
 174
 175
 176 void fastq_free(fastq_t* f)
 177 {
 178     free(f->buf);
 179     free(f);
 180 }
 181
 182
 183 typedef enum {
 184     FASTQ_STATE_ID1,  /* Reading ID1. */
 185     FASTQ_STATE_SEQ,  /* Reading the sequence. */
 186     FASTQ_STATE_ID2,  /* Reading ID2. */
 187     FASTQ_STATE_QUAL, /* Reading quality scores. */
 188 } fastq_parser_state_t;
 189
 190
 191 bool fastq_read(fastq_t* f, seq_t* seq)
 192 {
 193     seq->id1.n = seq->seq.n = seq->id2.n = seq->qual.n = 0;
 194     fastq_parser_state_t state = FASTQ_STATE_ID1;
 195     char* end = f->buf + f->readlen;
 196     do {
 197         while (f->next < end) {
 198             /* Consume pointless special characters prefixing IDs */
 199             if ((state == FASTQ_STATE_ID1 && f->linestart && f->next[0] == '@') ||
 200                 (state == FASTQ_STATE_ID2 && f->linestart && f->next[0] == '+')) {
 201                 f->linestart = false;
 202                 ++f->next;
 203                 continue;
 204             }
 205
 206             char* u = memchr(f->next, '\n', end - f->next);
 207             if (u == NULL) {
 208                 f->linestart = false;
 209                 u = end;
 210             }
 211             else f->linestart = true;
 212
 213             switch (state) {
 214                 case FASTQ_STATE_ID1:
 215                     str_append(&seq->id1, f->next, u - f->next);
 216                     if (f->linestart) state = FASTQ_STATE_SEQ;
 217                     break;
 218
 219                 case FASTQ_STATE_SEQ:
 220                     str_append(&seq->seq, f->next, u - f->next);
 221                     if (f->linestart) state = FASTQ_STATE_ID2;
 222                     break;
 223
 224                 case FASTQ_STATE_ID2:
 225                     str_append(&seq->id2, f->next, u - f->next);
 226                     if (f->linestart) state = FASTQ_STATE_QUAL;
 227                     break;
 228
 229                 case FASTQ_STATE_QUAL:
 230                     str_append(&seq->qual, f->next, u - f->next);
 231                     if (f->linestart) {
 232                         f->next = u + 1;
 233                         return true;
 234                     }
 235                     break;
 236             }
 237
 238             f->next = u + 1;
 239         }
 240
 241         /* Try to read more. */
 242         f->readlen = fread(f->buf, 1, parser_buf_size, f->file);
 243         f->next = f->buf;
 244         end = f->buf + f->readlen;
 245     } while (f->readlen);
 246
 247     return false;
 248 }
 249
 250
 251 void fastq_rewind(fastq_t* f)
 252 {
 253     rewind(f->file);
 254     f->next = f->buf;
 255     f->readlen = 0;
 256 }
 257
 258
 259 void fastq_print(FILE* fout, const seq_t* seq)
 260 {
 261     fprintf(fout, "@%s\n%s\n+%s\n%s\n",
 262                   seq->id1.s,
 263                   seq->seq.s,
 264                   seq->id2.s,
 265                   seq->qual.s );
 266 }
 267
 268