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);
  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(uint32_t seed, 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 = seed;
  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 static uint32_t seq_hash_seed = 0xc062fb4a;
 146
 147
 148 void seq_hash_set_seed(uint32_t seed)
 149 {
 150     seq_hash_seed = seed;
 151 }
 152
 153
 154 uint32_t seq_hash(const seq_t* seq)
 155 {
 156     uint32_t h = seq_hash_seed;
 157     h = murmurhash3(h, (uint8_t*) seq->id1.s, seq->id1.n);
 158     h = murmurhash3(h, (uint8_t*) seq->seq.s, seq->seq.n);
 159     h = murmurhash3(h, (uint8_t*) seq->id2.s, seq->id2.n);
 160     h = murmurhash3(h, (uint8_t*) seq->qual.s, seq->qual.n);
 161     return h;
 162 }
 163
 164
 165 static const size_t parser_buf_size = 1000000;
 166
 167
 168 struct fastq_t_
 169 {
 170     FILE* file;
 171     size_t readlen;
 172     char* buf;
 173     char* next;
 174     bool linestart;
 175 };
 176
 177
 178
 179 fastq_t* fastq_create(FILE* file)
 180 {
 181     fastq_t* f = malloc_or_die(sizeof(fastq_t));
 182     f->file = file;
 183     f->next = f->buf = malloc_or_die(parser_buf_size);
 184     f->readlen = 0;
 185     f->linestart = true;
 186     return f;
 187 }
 188
 189
 190 void fastq_free(fastq_t* f)
 191 {
 192     free(f->buf);
 193     free(f);
 194 }
 195
 196
 197 typedef enum {
 198     FASTQ_STATE_ID1,  /* Reading ID1. */
 199     FASTQ_STATE_SEQ,  /* Reading the sequence. */
 200     FASTQ_STATE_ID2,  /* Reading ID2. */
 201     FASTQ_STATE_QUAL, /* Reading quality scores. */
 202 } fastq_parser_state_t;
 203
 204
 205 bool fastq_read(fastq_t* f, seq_t* seq)
 206 {
 207     seq->id1.n = seq->seq.n = seq->id2.n = seq->qual.n = 0;
 208     fastq_parser_state_t state = FASTQ_STATE_ID1;
 209     char* end = f->buf + f->readlen;
 210     do {
 211         while (f->next < end) {
 212             /* Consume pointless special characters prefixing IDs */
 213             if ((state == FASTQ_STATE_ID1 && f->linestart && f->next[0] == '@') ||
 214                 (state == FASTQ_STATE_ID2 && f->linestart && f->next[0] == '+')) {
 215                 f->linestart = false;
 216                 ++f->next;
 217                 continue;
 218             }
 219
 220             char* u = memchr(f->next, '\n', end - f->next);
 221             if (u == NULL) {
 222                 f->linestart = false;
 223                 u = end;
 224             }
 225             else f->linestart = true;
 226
 227             switch (state) {
 228                 case FASTQ_STATE_ID1:
 229                     str_append(&seq->id1, f->next, u - f->next);
 230                     if (f->linestart) state = FASTQ_STATE_SEQ;
 231                     break;
 232
 233                 case FASTQ_STATE_SEQ:
 234                     str_append(&seq->seq, f->next, u - f->next);
 235                     if (f->linestart) state = FASTQ_STATE_ID2;
 236                     break;
 237
 238                 case FASTQ_STATE_ID2:
 239                     str_append(&seq->id2, f->next, u - f->next);
 240                     if (f->linestart) state = FASTQ_STATE_QUAL;
 241                     break;
 242
 243                 case FASTQ_STATE_QUAL:
 244                     str_append(&seq->qual, f->next, u - f->next);
 245                     if (f->linestart) {
 246                         f->next = u + 1;
 247                         return true;
 248                     }
 249                     break;
 250             }
 251
 252             f->next = u + 1;
 253         }
 254
 255         /* Try to read more. */
 256         f->readlen = fread(f->buf, 1, parser_buf_size, f->file);
 257         f->next = f->buf;
 258         end = f->buf + f->readlen;
 259     } while (f->readlen);
 260
 261     return false;
 262 }
 263
 264
 265 void fastq_rewind(fastq_t* f)
 266 {
 267     rewind(f->file);
 268     f->next = f->buf;
 269     f->readlen = 0;
 270 }
 271
 272
 273 void fastq_print(FILE* fout, const seq_t* seq)
 274 {
 275     fprintf(fout, "@%s\n%s\n+%s\n%s\n",
 276                   seq->id1.s,
 277                   seq->seq.s,
 278                   seq->id2.s,
 279                   seq->qual.s );
 280 }
 281
 282