src/fastq-kmers.c

   1
   2 /*
   3  * This file is part of fastq-tools.
   4  *
   5  * Copyright (c) 2011 by Daniel C. Jones <dcjones@cs.washington.edu>
   6  *
   7  * fastq-kmers :
   8  * Tabulate k-mer frequencies with FASTQ files.
   9  *
  10  */
  11
  12 #include "common.h"
  13 #include "parse.h"
  14 #include <stdlib.h>
  15 #include <stdio.h>
  16 #include <string.h>
  17 #include <stdint.h>
  18 #include <getopt.h>
  19 #include <zlib.h>
  20
  21
  22 #if defined(MSDOS) || defined(OS2) || defined(WIN32) || defined(__CYGWIN__)
  23 #  include <fcntl.h>
  24 #  include <io.h>
  25 #  define SET_BINARY_MODE(file) setmode(fileno(file), O_BINARY)
  26 #else
  27 #  define SET_BINARY_MODE(file)
  28 #endif
  29
  30 static const char* prog_name = "fastq-kmers";
  31
  32 void print_help()
  33 {
  34     fprintf( stderr,
  35 "fastq-kmers [OPTION]... [FILE]...\n"
  36 "Print kmer counts for the given kmer size.\n"
  37 "Output is in two tab-seperated columns for kmer and frequency.\n\n"
  38 "Options:\n"
  39 "  -k NUM, --size=NUM      kmer size (default: 1)\n"
  40 "  -h, --help              print this message\n"
  41 "  -V, --version           output version information and exit\n"
  42     );
  43 }
  44
  45 static int k;
  46
  47 int packkmer( const char* s, uint32_t* kmer, int k )
  48 {
  49     *kmer = 0;
  50     uint32_t nt;
  51     while (k--) {
  52         switch (s[k]) {
  53             case 'a':
  54             case 'A':
  55                 nt = 0;
  56                 break;
  57             case 'c':
  58             case 'C':
  59                 nt = 1;
  60                 break;
  61             case 'g':
  62             case 'G':
  63                 nt = 2;
  64                 break;
  65             case 't':
  66             case 'T':
  67                 nt = 3;
  68                 break;
  69
  70             default:
  71                 return 0;
  72         }
  73
  74         *kmer = (*kmer << 2) | nt;
  75     }
  76
  77     return 1;
  78 }
  79
  80
  81 void unpackkmer( uint32_t kmer, char* s, int k )
  82 {
  83     int i;
  84     uint32_t nt;
  85     for (i = 0; i < k; i++, s++) {
  86         nt = kmer & 0x3;
  87         kmer = kmer >> 2;
  88
  89         switch (nt) {
  90             case 0:
  91                 *s = 'A';
  92                 break;
  93             case 1:
  94                 *s = 'C';
  95                 break;
  96             case 2:
  97                 *s = 'G';
  98                 break;
  99             case 3:
 100                 *s = 'T';
 101                 break;
 102             default: break;
 103         }
 104     }
 105
 106     *s = '\0';
 107 }
 108
 109
 110 void count_fastq_kmers(FILE* fin, uint32_t* cs)
 111 {
 112     seq_t* seq = seq_create();
 113     fastq_t* fqf = fastq_create(fin);
 114     int i;
 115     int n;
 116     uint32_t kmer;
 117
 118     while (fastq_read(fqf, seq)) {
 119         n = (int)seq->seq.n - k + 1;
 120         for (i = 0; i < n; i++) {
 121             if( packkmer(seq->seq.s + i, &kmer, k) ) {
 122                 cs[kmer]++;
 123             }
 124         }
 125     }
 126
 127     seq_free(seq);
 128     fastq_free(fqf);
 129 }
 130
 131
 132 void print_kmer_freqs(FILE* fout, uint32_t* cs)
 133 {
 134     uint32_t n = 1 << (2*k); /* 4^k */
 135
 136     char* kmer_str = malloc((k+1)*sizeof(char));
 137     uint32_t kmer;
 138
 139     fprintf(fout, "kmer\tfrequency\n");
 140     for (kmer = 0; kmer < n; kmer++) {
 141         unpackkmer(kmer, kmer_str, k);
 142         fprintf(fout, "%s\t%u\n", kmer_str, cs[kmer]);
 143     }
 144
 145     free(kmer_str);
 146 }
 147
 148
 149 int main(int argc, char* argv[])
 150 {
 151     SET_BINARY_MODE(stdin);
 152     SET_BINARY_MODE(stdout);
 153
 154     k = 1;
 155
 156     uint32_t n;   /* number of kmers: 4^k */
 157     uint32_t* cs; /* counts */
 158
 159     FILE* fin;
 160
 161     int opt;
 162     int opt_idx;
 163     static struct option long_options[] =
 164         {
 165           {"size", no_argument,    0, 0},
 166           {"help", no_argument,    0, 'h'},
 167           {"version", no_argument, 0, 'V'},
 168           {0, 0, 0, 0}
 169         };
 170
 171     while (1) {
 172         opt = getopt_long(argc, argv, "k:hV", long_options, &opt_idx);
 173
 174         if( opt == -1 ) break;
 175
 176         switch (opt) {
 177             case 0:
 178                 if (long_options[opt_idx].flag != 0) break;
 179                 if (optarg) {
 180                     if( opt_idx == 1 ) {
 181                         k = atoi(optarg);
 182                     }
 183                 }
 184                 break;
 185
 186             case 'k':
 187                 k = atoi(optarg);
 188                 break;
 189
 190             case 'h':
 191                 print_help();
 192                 return 0;
 193
 194             case 'V':
 195                 print_version(stdout, prog_name);
 196                 return 0;
 197
 198             case '?':
 199                 return 1;
 200
 201             default:
 202                 abort();
 203         }
 204     }
 205
 206     if (k < 1) {
 207         fprintf(stderr, "Kmer size must be at least 1.");
 208         return 1;
 209     }
 210
 211     if (k > 16) {
 212         fprintf(stderr, "Kmer size must be at most 16.");
 213         return 1;
 214     }
 215
 216
 217     n = 1 << (2*k); /* i.e. 4^k */
 218     cs = malloc(n * sizeof(uint32_t));
 219     memset(cs, 0, n * sizeof(uint32_t));
 220
 221     if (cs == NULL) {
 222         fprintf(stderr, "Insufficient memory to tally kmers of size %d\n", k );
 223         return 1;
 224     }
 225
 226     if (optind >= argc || (argc - optind == 1 && strcmp(argv[optind],"-") == 0)) {
 227         count_fastq_kmers(stdin, cs);
 228     }
 229     else {
 230         for (; optind < argc; optind++) {
 231             fin = fopen(argv[optind], "rb");
 232             if (fin == NULL) {
 233                 fprintf(stderr, "No such file '%s'.\n", argv[optind]);
 234                 continue;
 235             }
 236
 237             count_fastq_kmers(fin, cs);
 238         }
 239     }
 240
 241     print_kmer_freqs( stdout, cs );
 242     free(cs);
 243
 244     return 0;
 245 }
 246
 247
 248
 249
 250
 251