src/fastq-qscale.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-qualscale:
   7  * Determine the quality score encoding used by a fastq file.
   8  *
   9  */
  10
  11 #include <getopt.h>
  12
  13 #include "common.h"
  14 #include "parse.h"
  15
  16 /* The scales here are taken from the wikipedia article for FASTQ, the relavent
  17  * part is reproduced here:
  18  *
  19  *  SSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSS.....................................................
  20  *  ..........................XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX......................
  21  *  ...............................IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII......................
  22  *  .................................JJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJ......................
  23  *  LLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLLL....................................................
  24  *  !"#$%&'()*+,-./0123456789:;<=>?@ABCDEFGHIJKLMNOPQRSTUVWXYZ[\]^_`abcdefghijklmnopqrstuvwxyz{|}~
  25  *  |                         |    |        |                              |                     |
  26  * 33                        59   64       73                            104                   126
  27  *
  28  *    S - Sanger        Phred+33,  raw reads typically (0, 40)
  29  *    X - Solexa        Solexa+64, raw reads typically (-5, 40)
  30  *    I - Illumina 1.3+ Phred+64,  raw reads typically (0, 40)
  31  *    J - Illumina 1.5+ Phred+64,  raw reads typically (3, 40)
  32  *       with 0=unused, 1=unused, 2=Read Segment Quality Control Indicator (bold)
  33  *       (Note: See discussion above).
  34  *        L - Illumina 1.8+ Phred+33,  raw reads typically (0, 41)
  35  */
  36
  37
  38 typedef struct qual_scale_t_
  39 {
  40     const char* description;
  41     char min_qual, max_qual;
  42 } qual_scale_t;
  43
  44
  45 /* When the scale is ambiguous, we choose the first compatible one. Hence these
  46  * are ordered roughly by increasing exclusivity. */
  47 #define NUM_SCALES 5
  48 static qual_scale_t scales[NUM_SCALES] =
  49 {
  50     {"Sanger/Phred+33",       '!', 'I'},
  51     {"Illumina 1.8/Phred+33", '!', 'J'},
  52     {"Illumina 1.5/Phred+64", 'B', 'h'},
  53     {"Illumina 1.3/Phred+64", '@', 'h'},
  54     {"Solexa/Solexa+64",      ';', 'h'}
  55 };
  56
  57
  58 /* Return true if x has excatly one 1 bit. */
  59 bool single_bit(uint32_t x)
  60 {
  61     return x && !(x & (x - 1));
  62 }
  63
  64
  65 /* Make a bitset of compatible scales. */
  66 uint32_t make_bitset(char min_qual, char max_qual)
  67 {
  68     uint32_t s = 0;
  69     uint32_t i;
  70     for (i = 0; i < NUM_SCALES; ++i) {
  71         if (scales[i].min_qual <= min_qual &&
  72             scales[i].max_qual >= max_qual) {
  73             s |= 1 << i;
  74         }
  75     }
  76
  77     return s;
  78 }
  79
  80
  81 void fastq_qualscale(const char* fn, FILE* fin)
  82 {
  83     char min_qual = '~', max_qual = '!';
  84
  85     fastq_t* fqf = fastq_create(fin);
  86     seq_t* seq = seq_create();
  87
  88     /* Scales compatible with the data so far. */
  89     uint32_t compat_scales;
  90
  91     size_t n = 100000;
  92
  93     while (n-- && fastq_read(fqf, seq) && !single_bit(compat_scales)) {
  94         char* q = seq->qual.s;
  95         while (*q) {
  96             if (*q < min_qual) min_qual = *q;
  97             if (*q > max_qual) max_qual = *q;
  98             ++q;
  99         }
 100
 101         compat_scales = make_bitset(min_qual, max_qual);
 102         if (compat_scales == 0 || single_bit(compat_scales)) break;
 103     }
 104
 105     seq_free(seq);
 106     fastq_free(fqf);
 107
 108     if (compat_scales == 0) {
 109         printf("%s: Unknown scale ['%c', '%c']\n", fn, min_qual, max_qual);
 110     }
 111     else {
 112         /* low order bit */
 113         unsigned int i;
 114         for (i = 0; !(compat_scales & (1 << i)); ++i) {}
 115         printf("%s: %s\n", fn, scales[i].description);
 116     }
 117 }
 118
 119
 120 static const char* prog_name = "fastq-qscale";
 121
 122
 123 void print_help()
 124 {
 125     fprintf(stderr,
 126 "fastq-qscale [OPTION] [FILE]...\n"
 127 "Detect and output the quality score scale for each file given as an argument.\n\n"
 128 "Options:\n"
 129 "  -h, --help              print this message\n"
 130 "  -V, --version           output version information and exit\n"
 131     );
 132 }
 133
 134
 135 int main(int argc, char* argv[])
 136 {
 137     static struct option long_options[] =
 138     {
 139         {"help",    no_argument, NULL, 'h'},
 140         {"version", no_argument, NULL, 'V'},
 141         {0, 0, 0, 0}
 142     };
 143
 144     int opt, opt_idx;
 145     while (true) {
 146         opt = getopt_long(argc, argv, "hV", long_options, &opt_idx);
 147         if (opt == -1) break;
 148
 149         switch (opt) {
 150             case 'h':
 151                 print_help();
 152                 return 0;
 153
 154             case 'V':
 155                 print_version(stdout, prog_name);
 156                 return 0;
 157
 158             case '?':
 159                 return 1;
 160
 161             default:
 162                 abort();
 163         }
 164     }
 165
 166     if (optind >= argc) {
 167         fastq_qualscale("stdin", stdin);
 168     }
 169     else {
 170         for (; optind < argc; optind++) {
 171             FILE* fin = fopen(argv[optind], "rb");
 172             if (fin == NULL) {
 173                 fprintf(stderr, "No such file '%s'.\n", argv[optind]);
 174                 continue;
 175             }
 176
 177             fastq_qualscale(argv[optind], fin);
 178             fclose(fin);
 179         }
 180     }
 181
 182     return EXIT_SUCCESS;
 183 }
 184