src/fastq-grep.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-grep :
   8  * Regular expression searches of the sequences within a FASTQ file.
   9  *
  10  */
  11
  12
  13 #include "common.h"
  14 #include "parse.h"
  15 #include <stdio.h>
  16 #include <string.h>
  17 #include <getopt.h>
  18 #include <zlib.h>
  19 #include <pcre.h>
  20
  21
  22 #if defined(MSDOS) || defined(OS2) || defined(WIN32) || defined(__CYGWIN__)
  23 #  include <fcntl.h>
  24 #  include <io.h>
  25 #  define SET_BINARY_MODE(file) setmode(fileno(file), O_BINARY)
  26 #else
  27 #  define SET_BINARY_MODE(file)
  28 #endif
  29
  30
  31 static const char* prog_name = "fastq-grep";
  32
  33
  34 void print_help()
  35 {
  36     fprintf(stdout,
  37 "fastq-grep [OPTION]... PATTERN [FILE]...\n"
  38 "Search for PATTERN in the read sequences in each FILE or standard input.\n"
  39 "PATTERN, by default, is a perl compatible regular expression.\n\n"
  40 "Options:\n"
  41 "  -i, --id                match the read id (by default, sequence is matched)\n"
  42 "  -v, --invert-match      select nonmatching entries\n"
  43 "  -m, --mismatches=FILE   output mismatching entries to the given file\n"
  44 "  -c, --count             output only the number of matching sequences\n"
  45 "  -h, --help              print this message\n"
  46 "  -V, --version           output version information and exit\n"
  47     );
  48 }
  49
  50 static int invert_flag;
  51 static int count_flag;
  52 static int id_flag;
  53
  54
  55
  56 void fastq_grep(FILE* fin, FILE* fout, FILE* mismatch_file, pcre* re)
  57 {
  58     int rc;
  59     int ovector[3];
  60     size_t count = 0;
  61
  62     fastq_t* fqf = fastq_create(fin);
  63     seq_t* seq = seq_create();
  64
  65     while (fastq_read(fqf, seq)) {
  66
  67         rc = pcre_exec(re,          /* pattern */
  68                        NULL,        /* extra data */
  69                        id_flag ? seq->id1.s : seq->seq.s,
  70                        id_flag ? seq->id1.n : seq->seq.n,
  71                        0,           /* subject offset */
  72                        0,           /* options */
  73                        ovector,     /* output vector */
  74                        3         ); /* output vector length */
  75
  76         if ((invert_flag && rc == PCRE_ERROR_NOMATCH) || (!invert_flag && rc >= 0)) {
  77             if (count_flag) count++;
  78             else            fastq_print(fout, seq);
  79         }
  80         else if (mismatch_file) {
  81             fastq_print(mismatch_file, seq);
  82         }
  83     }
  84
  85     seq_free(seq);
  86     fastq_free(fqf);
  87
  88     if (count_flag) fprintf(fout, "%zu\n", count);
  89 }
  90
  91
  92
  93 int main(int argc, char* argv[])
  94 {
  95     SET_BINARY_MODE(stdin);
  96     SET_BINARY_MODE(stdout);
  97
  98     const char* pat;
  99     pcre* re;
 100     const char* pat_error;
 101     int pat_error_offset;
 102
 103     FILE*  fin;
 104
 105
 106     invert_flag  = 0;
 107     count_flag   = 0;
 108     id_flag      = 0;
 109
 110     int opt;
 111     int opt_idx;
 112
 113     FILE* mismatch_file = NULL;
 114
 115     static struct option long_options[] =
 116         {
 117           {"id",           no_argument, &id_flag,     1},
 118           {"invert-match", no_argument, &invert_flag, 1},
 119           {"mismatches",   required_argument, NULL, 'm'},
 120           {"count",        no_argument, &count_flag,  1},
 121           {"help",         no_argument, NULL, 'h'},
 122           {"version",      no_argument, NULL, 'V'},
 123           {0, 0, 0, 0}
 124         };
 125
 126     while (1) {
 127         opt = getopt_long(argc, argv, "ivmchV", long_options, &opt_idx);
 128
 129         if (opt == -1) break;
 130
 131         switch (opt) {
 132             case 0:
 133                 if (long_options[opt_idx].flag != 0) break;
 134                 if (optarg) {
 135                 }
 136                 break;
 137
 138             case 'i':
 139                 id_flag = 1;
 140                 break;
 141
 142             case 'v':
 143                 invert_flag = 1;
 144                 break;
 145
 146             case 'm':
 147                 mismatch_file = fopen(optarg, "w");
 148                 if (mismatch_file == NULL) {
 149                     fprintf(stderr, "No such file '%s'.\n", optarg);
 150                     return 1;
 151                 }
 152                 break;
 153
 154             case 'c':
 155                 count_flag = 1;
 156                 break;
 157
 158             case 'h':
 159                 print_help();
 160                 return 0;
 161
 162             case 'V':
 163                 print_version(stdout, prog_name);
 164                 return 0;
 165
 166             case '?':
 167                 return 1;
 168
 169             default:
 170                 abort();
 171         }
 172     }
 173
 174     if (optind >= argc) {
 175         fprintf(stderr, "A pattern must be specified.\n");
 176         return 1;
 177     }
 178
 179     pat = argv[optind++];
 180     re = pcre_compile( pat, PCRE_CASELESS, &pat_error, &pat_error_offset, NULL );
 181
 182
 183     if (re == NULL) {
 184         fprintf(stderr, "Syntax error in PCRE pattern at offset: %d: %s\n",
 185                 pat_error_offset, pat_error );
 186         return 1;
 187     }
 188
 189
 190     if (optind >= argc || (argc - optind == 1 && strcmp(argv[optind],"-") == 0)) {
 191         fastq_grep(stdin, stdout, mismatch_file, re);
 192     }
 193     else {
 194         for (; optind < argc; optind++) {
 195             fin = fopen(argv[optind], "rb");
 196             if (fin == NULL) {
 197                 fprintf(stderr, "No such file '%s'.\n", argv[optind]);
 198                 continue;
 199             }
 200
 201             fastq_grep(fin, stdout, mismatch_file, re);
 202
 203             fclose(fin);
 204         }
 205     }
 206
 207     pcre_free(re);
 208     if (mismatch_file) fclose(mismatch_file);
 209
 210     return 0;
 211 }
 212
 213
 214