src/fastq-kmers.c

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