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