src/fastq-uniq.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-uniq :
   8  * Collapsing a fastq file into only unique read sequences.
   9  */
  10
  11
  12 #include "common.h"
  13 #include "hash.h"
  14 #include "parse.h"
  15 #include <string.h>
  16 #include <zlib.h>
  17 #include <getopt.h>
  18
  19
  20 #if defined(MSDOS) || defined(OS2) || defined(WIN32) || defined(__CYGWIN__)
  21 #  include <fcntl.h>
  22 #  include <io.h>
  23 #  define SET_BINARY_MODE(file) setmode(fileno(file), O_BINARY)
  24 #else
  25 #  define SET_BINARY_MODE(file)
  26 #endif
  27
  28
  29 static int help_flag;
  30 static int verbose_flag;
  31 size_t total_reads;
  32
  33 void print_help()
  34 {
  35     fprintf(stderr,
  36 "fastq-uniq [OPTION] [FILE]...\n"
  37 "Output a non-redundant FASTQ file, in which there are no duplicate reads.\n"
  38 "(Warning: this program can be somewhat memory intensive.)\n\n"
  39 "Options:\n"
  40 "  -h, --help              print this message\n"
  41 "  -v, --verbose           print status along the way\n"
  42     );
  43 }
  44
  45
  46 void fastq_hash(FILE* fin, hash_table* T)
  47 {
  48     fastq_t* fqf = fastq_open(fin);
  49     seq_t* seq = fastq_alloc_seq();
  50
  51     while (fastq_next(fqf, seq)) {
  52         inc_hash_table(T, seq->seq.s, seq->seq.n);
  53
  54         total_reads++;
  55         if (verbose_flag && total_reads % 100000 == 0) {
  56             fprintf(stderr, "%zu reads processed ...\n", total_reads);
  57         }
  58     }
  59
  60     fastq_free_seq(seq);
  61     fastq_close(fqf);
  62 }
  63
  64
  65 int compare_hashed_value_count(const void* x, const void* y)
  66 {
  67     hashed_value* const * a = x;
  68     hashed_value* const * b = y;
  69
  70     if( (*a)->count > (*b)->count ) return -1;
  71     if( (*a)->count < (*b)->count ) return 1;
  72     return 0;
  73 }
  74
  75
  76
  77 void print_hash_table(FILE* fout, hash_table* T)
  78 {
  79     hashed_value** S = dump_hash_table(T);
  80     qsort(S, T->m, sizeof(hashed_value*), compare_hashed_value_count);
  81
  82     size_t i;
  83     for (i = 0; i < T->m; i++) {
  84         fprintf(fout, ">unique-read-%07zu (%zu copies)\n", i, S[i]->count);
  85         fwrite(S[i]->value, S[i]->len, sizeof(char), fout);
  86         fprintf(fout, "\n");
  87     }
  88     free(S);
  89 }
  90
  91
  92
  93 int main(int argc, char* argv[])
  94 {
  95     SET_BINARY_MODE(stdin);
  96     SET_BINARY_MODE(stdout);
  97
  98     hash_table* T = create_hash_table();
  99
 100     FILE* fin   ;
 101
 102     help_flag = 0;
 103
 104     int opt;
 105     int opt_idx;
 106
 107     static struct option long_options[] =
 108     {
 109         {"help",    no_argument, &help_flag,    1},
 110         {"verbose", no_argument, &verbose_flag, 1},
 111         {0, 0, 0, 0}
 112     };
 113
 114     while (1) {
 115         opt = getopt_long(argc, argv, "hv", long_options, &opt_idx);
 116
 117         if (opt == -1) break;
 118
 119         switch (opt) {
 120             case 0:
 121                 if (long_options[opt_idx].flag != 0) break;
 122                 if (optarg) {
 123                 }
 124                 break;
 125
 126             case 'h':
 127                 help_flag = 1;
 128                 break;
 129
 130             case 'v':
 131                 verbose_flag = 1;
 132                 break;
 133
 134             case '?':
 135                 return 1;
 136
 137             default:
 138                 abort();
 139         }
 140     }
 141
 142     if (help_flag) {
 143         print_help();
 144         return 0;
 145     }
 146
 147     if (optind >= argc || (argc - optind == 1 && strcmp(argv[optind],"-") == 0)) {
 148         fastq_hash(stdin, T);
 149     }
 150     else {
 151         for (; optind < argc; optind++) {
 152             fin = fopen(argv[optind], "rb");
 153             if (fin == NULL) {
 154                 fprintf(stderr, "No such file '%s'.\n", argv[optind]);
 155                 continue;
 156             }
 157
 158             fastq_hash(fin, T);
 159         }
 160     }
 161
 162     print_hash_table(stdout, T);
 163
 164     destroy_hash_table(T);
 165     return 0;
 166 }