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 "fastq-common.h"
  13 #include "fastq-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
  31 void print_help()
  32 {
  33     fprintf( stderr,
  34 "fastq-kmers [OPTION]... [FILE]...\n"
  35 "Print kmer counts for the given kmer size.\n"
  36 "Output is in two tab-seperated columns for kmer and frequency.\n\n"
  37 "Options:\n"
  38 "  -h, --help              print this message\n"
  39 "  -k, --size              kmer size (default: 1)\n"
  40     );
  41 }
  42
  43 static int help_flag;
  44 static int k;
  45
  46 int packkmer( const char* s, uint32_t* kmer, int k )
  47 {
  48     *kmer = 0;
  49     uint32_t nt;
  50     while (k--) {
  51         switch (s[k]) {
  52             case 'a':
  53             case 'A':
  54                 nt = 0;
  55                 break;
  56             case 'c':
  57             case 'C':
  58                 nt = 1;
  59                 break;
  60             case 'g':
  61             case 'G':
  62                 nt = 2;
  63                 break;
  64             case 't':
  65             case 'T':
  66                 nt = 3;
  67                 break;
  68
  69             default:
  70                 return 0;
  71         }
  72
  73         *kmer = (*kmer << 2) | nt;
  74     }
  75
  76     return 1;
  77 }
  78
  79
  80 void unpackkmer( uint32_t kmer, char* s, int k )
  81 {
  82     int i;
  83     uint32_t nt;
  84     for (i = 0; i < k; i++, s++) {
  85         nt = kmer & 0x3;
  86         kmer = kmer >> 2;
  87
  88         switch (nt) {
  89             case 0:
  90                 *s = 'A';
  91                 break;
  92             case 1:
  93                 *s = 'C';
  94                 break;
  95             case 2:
  96                 *s = 'G';
  97                 break;
  98             case 3:
  99                 *s = 'T';
 100                 break;
 101             default: break;
 102         }
 103     }
 104
 105     *s = '\0';
 106 }
 107
 108
 109 void count_fastq_kmers(FILE* fin, uint32_t* cs)
 110 {
 111     seq_t* seq = fastq_alloc_seq();
 112     fastq_t* fqf = fastq_open(fin);
 113     int i;
 114     int n;
 115     uint32_t kmer;
 116
 117     while (fastq_next(fqf, seq)) {
 118         n = (int)seq->seq.n - k + 1;
 119         for (i = 0; i < n; i++) {
 120             if( packkmer(seq->seq.s + i, &kmer, k) ) {
 121                 cs[kmer]++;
 122             }
 123         }
 124     }
 125
 126     fastq_free_seq(seq);
 127     fastq_close(fqf);
 128 }
 129
 130
 131 void print_kmer_freqs(FILE* fout, uint32_t* cs)
 132 {
 133     uint32_t n = 1 << (2*k); /* 4^k */
 134
 135     char* kmer_str = malloc((k+1)*sizeof(char));
 136     uint32_t kmer;
 137
 138     fprintf(fout, "kmer\tfrequency\n");
 139     for (kmer = 0; kmer < n; kmer++) {
 140         unpackkmer(kmer, kmer_str, k);
 141         fprintf(fout, "%s\t%u\n", kmer_str, cs[kmer]);
 142     }
 143
 144     free(kmer_str);
 145 }
 146
 147
 148 int main(int argc, char* argv[])
 149 {
 150     SET_BINARY_MODE(stdin);
 151     SET_BINARY_MODE(stdout);
 152
 153     help_flag = 0;
 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           {"help", no_argument, &help_flag, 1},
 166           {"size", no_argument, 0, 0},
 167           {0, 0, 0, 0}
 168         };
 169
 170     while (1) {
 171         opt = getopt_long(argc, argv, "hk:", long_options, &opt_idx);
 172
 173         if( opt == -1 ) break;
 174
 175         switch (opt) {
 176             case 0:
 177                 if (long_options[opt_idx].flag != 0) break;
 178                 if (optarg) {
 179                     if( opt_idx == 1 ) {
 180                         k = atoi(optarg);
 181                     }
 182                 }
 183                 break;
 184
 185             case 'h':
 186                 help_flag = 1;
 187                 break;
 188
 189             case 'k':
 190                 k = atoi(optarg);
 191                 break;
 192
 193             case '?':
 194                 return 1;
 195
 196             default:
 197                 abort();
 198         }
 199     }
 200
 201     if (help_flag) {
 202         print_help();
 203         return 0;
 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