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