src/fastq-sort.c

   1 /*
   2  * This file is part of fastq-tools.
   3  *
   4  * Copyright (c) 2012 by Daniel C. Jones <dcjones@cs.washington.edu>
   5  *
   6  * fastq-sample :
   7  * Sample reads with or without replacement from a FASTQ file.
   8  *
   9  */
  10
  11 #include <assert.h>
  12 #include <getopt.h>
  13 #include <string.h>
  14
  15 #include "common.h"
  16 #include "parse.h"
  17
  18
  19 typedef struct seq_array_t_
  20 {
  21     seq_t* seqs;
  22
  23     /* Number of seq objects. */
  24     size_t n;
  25
  26     /* Size reserved in seqs. */
  27     size_t size;
  28
  29     /* Space reserved for strings. */
  30     char* data;
  31
  32     /* Data used. */
  33     size_t data_used;
  34
  35     /* Total data size. */
  36     size_t data_size;
  37 } seq_array_t;
  38
  39
  40 seq_array_t* seq_array_create(size_t data_size)
  41 {
  42     seq_array_t* a = malloc_or_die(sizeof(seq_array_t));
  43     a->size = 1024;
  44     a->n = 0;
  45     a->seqs = malloc_or_die(sizeof(seq_t));
  46
  47     a->data_size = data_size;
  48     a->data_used = 0;
  49     a->data = malloc_or_die(data_size);
  50
  51     return a;
  52 }
  53
  54
  55 void seq_array_free(seq_array_t* a)
  56 {
  57     free(a->seqs);
  58     free(a->data);
  59     free(a);
  60 }
  61
  62
  63 /* Push a fastq entry to back of the array. Return false if there was not enough
  64  * space. */
  65 bool seq_array_push(seq_array_t* a, const seq_t* seq)
  66 {
  67     size_t size_needed = (seq->id1.n + 1) + (seq->seq.n + 1) +
  68                          (seq->id2.n + 1) + (seq->qual.n + 1);
  69
  70     if (size_needed > a->data_size - a->data_used) return false;
  71
  72     if (a->n == a->size) {
  73         a->size *= 2;
  74         a->seqs = realloc_or_die(a->seqs, a->size * sizeof(seq_t));
  75     }
  76
  77     memcpy(seq->id1.s, &a->data[a->data_used], seq->id1.n + 1);
  78     a->seqs[a->n].id1.s = &a->data[a->data_used];
  79     a->seqs[a->n].id1.n = seq->id1.n + 1;
  80     a->data_used += seq->id1.n + 1;
  81
  82     memcpy(seq->seq.s, &a->data[a->data_used], seq->seq.n + 1);
  83     a->seqs[a->n].seq.s = &a->data[a->data_used];
  84     a->seqs[a->n].seq.n = seq->seq.n + 1;
  85     a->data_used += seq->seq.n + 1;
  86
  87     memcpy(seq->id2.s, &a->data[a->data_used], seq->id2.n + 1);
  88     a->seqs[a->n].id2.s = &a->data[a->data_used];
  89     a->seqs[a->n].id2.n = seq->id2.n + 1;
  90     a->data_used += seq->id2.n + 1;
  91
  92     memcpy(seq->qual.s, &a->data[a->data_used], seq->qual.n + 1);
  93     a->seqs[a->n].qual.s = &a->data[a->data_used];
  94     a->seqs[a->n].qual.n = seq->qual.n + 1;
  95     a->data_used += seq->qual.n + 1;
  96
  97     ++a->n;
  98
  99     return true;
 100 }
 101
 102
 103 void seq_array_clear(seq_array_t* a)
 104 {
 105     a->n = 0;
 106     a->data_used = 0;
 107 }
 108
 109
 110 void seq_array_sort(seq_array_t* a, int (*cmp)(const void*, const void*))
 111 {
 112     qsort(a->seqs, a->n, sizeof(seq_t), cmp);
 113 }
 114
 115
 116 int seq_cmp_hash(const void* a_, const void* b_)
 117 {
 118     const seq_t* a = (seq_t*) a_;
 119     const seq_t* b = (seq_t*) b_;
 120     /* TODO: hash and compare. */
 121     return 0;
 122 }
 123
 124
 125 static const char* prog_name = "fastq-sort";
 126
 127
 128 void print_help()
 129 {
 130     fprintf(stdout,
 131 "fastq-sort [OPTION]... [FILE]...\n"
 132 "Concatenate and sort FASTQ files and write to standard output.\n"
 133 "Options:\n"
 134 "  -h, --help      print this message\n"
 135 "  -V, --version   output version information and exit\n"
 136    );
 137 }
 138
 139
 140
 141 int main(int argc, char* argv[])
 142 {
 143     int opt, opt_idx;
 144     size_t buffer_size = 100000000;
 145     int (*cmp)(const void*, const void*) = seq_cmp_hash;;
 146
 147     static struct option long_options[] =
 148     {
 149         {"help",    no_argument, NULL, 'h'},
 150         {"version", no_argument, NULL, 'V'},
 151         {0, 0, 0, 0}
 152     };
 153
 154     while (true) {
 155         opt = getopt_long(argc, argv, "hV", long_options, &opt_idx);
 156         if (opt == -1) break;
 157
 158         switch (opt) {
 159             case 'h':
 160                 print_help();
 161                 return 0;
 162
 163             case 'V':
 164                 print_version(stdout, prog_name);
 165                 return 0;
 166
 167             case '?':
 168                 return 1;
 169
 170             default:
 171                 abort();
 172         }
 173     }
 174
 175     seq_array_t* a = seq_array_create(buffer_size);
 176     seq_t* seq = seq_create();
 177
 178     fastq_t* f;
 179     if (optind >= argc) {
 180         f = fastq_create(stdin);
 181         while (fastq_read(f, seq)) {
 182             if (!seq_array_push(a, seq)) {
 183                 seq_array_sort(a, cmp);
 184
 185                 /* TODO: dump a to a temporary file. Push that file name to an
 186                  * array somewhere.
 187                  * */
 188
 189                 seq_array_clear(a);
 190                 if (seq_array_push(a, seq)) {
 191                     fprintf(stderr, "The buffer size is to small.\n");
 192                     return EXIT_FAILURE;
 193                 }
 194             }
 195         }
 196         fastq_free(f);
 197     }
 198     else {
 199         FILE* file;
 200         for (; optind < argc; ++optind) {
 201             file = fopen(argv[optind], "rb");
 202             if (file == NULL) {
 203                 fprintf(stderr, "Cannot open %s for reading.\n", argv[optind]);
 204                 return EXIT_FAILURE;
 205             }
 206             f = fastq_create(file);
 207
 208             while (fastq_read(f, seq)) {
 209                 if (!seq_array_push(a, seq)) {
 210                     seq_array_sort(a, cmp);
 211
 212                     /* TODO: dump a to a temporary file. Push that file name to
 213                      * an array somewhere.  */
 214
 215                     seq_array_clear(a);
 216                     if (seq_array_push(a, seq)) {
 217                         fprintf(stderr, "The buffer size is to small.\n");
 218                         return EXIT_FAILURE;
 219                     }
 220                 }
 221             }
 222
 223             fastq_free(f);
 224             fclose(file);
 225         }
 226     }
 227
 228     if (a->n > 0) {
 229         seq_array_sort(a, cmp);
 230
 231         /* TODO: special case if everything fit into a. Just dump it to output.
 232          */
 233
 234         /* TODO: dump to a temp file, push file name to the stack. */
 235     }
 236
 237     /* TODO: Merge sort on the external files. */
 238
 239     seq_free(seq);
 240     seq_array_free(a);
 241
 242     return EXIT_SUCCESS;
 243 }
 244
 245