tests/random_fastq.c

   1
   2 /*
   3     random_fastq
   4     ------------
   5     Generate random data in FASTQ format.
   6 */
   7
   8 #include <stdlib.h>
   9 #include <stdio.h>
  10 #include <getopt.h>
  11
  12 static void print_help()
  13 {
  14     printf(
  15 "Usage: random_fastq [option]...\n"
  16 "Generate an endless stream of random FASTQ data to standard out.\n\n"
  17 "Options:\n"
  18 "    TODO\n\n"
  19 "Beware: the only purpose of this program is test quip.\n"
  20 "No particular guarantees are made.\n\n");
  21 }
  22
  23
  24 /* Draw n samples from a categorial distribution of size k with cumulative
  25  * distribution given by cs and element us. The sample is stored in xs which
  26  * is assumed to be of appropriate length. */
  27 void randcat(const char* us, const double* cs, size_t k, char* xs, size_t n)
  28 {
  29     size_t i;
  30     double r;
  31     for (i = 0; i < n; ++i) {
  32         r = drand48();
  33         int a = 0, b = (int) k, c;
  34         do {
  35             c = (a + b) / 2;
  36             if (r <= cs[c]) b = c;
  37             else            a = c + 1;
  38         } while (a < b);
  39
  40         xs[i] = us[a];
  41     }
  42 }
  43
  44
  45 int main(int argc, char* argv[])
  46 {
  47     static struct option long_options[] =
  48     {
  49         {"min-length", required_argument, NULL, 'm'},
  50         {"max-length", required_argument, NULL, 'M'},
  51         {"length",     required_argument, NULL, 'l'},
  52         {"id-length",  required_argument, NULL, 'i'},
  53         {"help",       no_argument,       NULL, 'h'},
  54         {NULL, 0, NULL, 0}
  55     };
  56
  57     size_t min_len = 100, max_len = 100;
  58     size_t id_len = 50;
  59
  60     int opt, opt_idx;
  61     while (1) {
  62         opt = getopt_long(argc, argv, "h", long_options, &opt_idx);
  63
  64         if (opt == -1) break;
  65
  66         switch (opt) {
  67             case 'm':
  68                 min_len = (size_t) strtoul(optarg, NULL, 10);
  69                 break;
  70
  71             case 'M':
  72                 max_len = (size_t) strtoul(optarg, NULL, 10);
  73                 break;
  74
  75             case 'l':
  76                 min_len = max_len = (size_t) strtoul(optarg, NULL, 10);
  77                 break;
  78
  79             case 'i':
  80                 id_len = (size_t) strtoul(optarg, NULL, 10);
  81                 break;
  82
  83             case 'h':
  84                 print_help();
  85                 return EXIT_SUCCESS;
  86
  87             case '?':
  88                 return EXIT_FAILURE;
  89
  90             default:
  91                 abort();
  92         }
  93     }
  94
  95     char nucleotides[5] = {'A', 'C', 'G', 'T', 'N'};
  96     double nuc_cs[5] = {0.28, 0.49, 0.70, 0.90, 1.00};
  97
  98     char qualities[64];
  99     double qual_cs[64];
 100     size_t i;
 101     double last_c = 0.0;
 102     for (i = 0; i < 64; ++i) {
 103         qualities[i] = '!' + i;
 104         last_c = qual_cs[i] = last_c + 1.0 / 64.0;
 105     }
 106
 107     char id_chars[94];
 108     double id_cs[94];
 109     last_c = 0.0;
 110     for (i = 0; i < 94; ++i) {
 111         id_chars[i] = '!' + i;
 112         last_c = id_cs[i] = last_c + 1.0 / 94.0;
 113     }
 114
 115     char* id   = malloc(id_len + 1);
 116     char* seq  = malloc(max_len + 1);
 117     char* qual = malloc(max_len + 1);
 118     size_t len = min_len;
 119
 120     while (1) {
 121         if (max_len > min_len) {
 122             len = min_len + (size_t) (drand48() * (double) (1 + max_len - min_len));
 123         }
 124
 125         randcat(id_chars, id_cs, 94, id, id_len);
 126         id[id_len] = '\0';
 127
 128         randcat(nucleotides, nuc_cs, 5, seq, len);
 129         seq[len] = '\0';
 130
 131         randcat(qualities, qual_cs, 64, qual, len);
 132         qual[len] = '\0';
 133
 134         printf(
 135             "@%s\n%s\n+\n%s\n",
 136             id, seq, qual);
 137     }
 138
 139     /* Yeah, right. As if we'll ever get here. */
 140     free(id);
 141     free(seq);
 142     free(qual);
 143
 144     return EXIT_SUCCESS;
 145 }
 146