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     help_flag = 0;
 148     k = 1;
 149
 150     uint32_t n;   /* number of kmers: 4^k */
 151     uint32_t* cs; /* counts */
 152
 153     FILE* fin;
 154     gzFile gzfin;
 155
 156     int opt;
 157     int opt_idx;
 158     static struct option long_options[] =
 159         {
 160           {"help", no_argument, &help_flag, 1},
 161           {"size", no_argument, 0, 0},
 162           {0, 0, 0, 0}
 163         };
 164
 165     while (1) {
 166         opt = getopt_long(argc, argv, "hk:", long_options, &opt_idx);
 167
 168         if( opt == -1 ) break;
 169
 170         switch (opt) {
 171             case 0:
 172                 if (long_options[opt_idx].flag != 0) break;
 173                 if (optarg) {
 174                     if( opt_idx == 1 ) {
 175                         k = atoi(optarg);
 176                     }
 177                 }
 178                 break;
 179
 180             case 'h':
 181                 help_flag = 1;
 182                 break;
 183
 184             case 'k':
 185                 k = atoi(optarg);
 186                 break;
 187
 188             case '?':
 189                 return 1;
 190
 191             default:
 192                 abort();
 193         }
 194     }
 195
 196     if (help_flag) {
 197         print_help();
 198         return 0;
 199     }
 200
 201     if (k < 1) {
 202         fprintf(stderr, "Kmer size must be at least 1.");
 203         return 1;
 204     }
 205
 206     if (k > 16) {
 207         fprintf(stderr, "Kmer size must be at most 16.");
 208         return 1;
 209     }
 210
 211
 212     n = 1 << (2*k); /* i.e. 4^k */
 213     cs = malloc(n * sizeof(uint32_t));
 214     memset(cs, 0, n * sizeof(uint32_t));
 215
 216     if (cs == NULL) {
 217         fprintf(stderr, "Insufficient memory to tally kmers of size %d\n", k );
 218         return 1;
 219     }
 220
 221     if (optind >= argc || (argc - optind == 1 && strcmp(argv[optind],"-") == 0)) {
 222         gzfin = gzdopen( fileno(stdin), "rb" );
 223         if (gzfin == NULL) {
 224             fprintf(stderr, "Malformed file 'stdin'.\n");
 225             return 1;
 226         }
 227
 228         count_fastq_kmers(gzfin, cs);
 229
 230         gzclose(gzfin);
 231     }
 232     else {
 233         for (; optind < argc; optind++) {
 234             fin = fopen(argv[optind], "rb");
 235             if (fin == NULL) {
 236                 fprintf(stderr, "No such file '%s'.\n", argv[optind]);
 237                 continue;
 238             }
 239
 240             gzfin = gzdopen(fileno(fin), "rb");
 241             if (gzfin == NULL) {
 242                 fprintf(stderr, "Malformed file '%s'.\n", argv[optind]);
 243                 continue;
 244             }
 245
 246             count_fastq_kmers(gzfin, cs);
 247
 248             gzclose(gzfin);
 249         }
 250     }
 251
 252     print_kmer_freqs( stdout, cs );
 253     free(cs);
 254
 255     return 0;
 256 }
 257
 258
 259
 260
 261
 262