src/fastq-grep.c

   1
   2
   3 /*
   4  * fastq-grep: regular expression searches of the sequences within a fastq file
   5  *
   6  * Febuary 2011 / Daniel Jones <dcjones@cs.washington.edu>
   7  *
   8  */
   9
  10
  11 #include <stdio.h>
  12 #include <string.h>
  13 #include <getopt.h>
  14 #include <zlib.h>
  15 #include <pcre.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-grep [OPTION]... PATTERN [FILE]...\n"
  34 "Search for PATTERN in the read sequences in each FILE or standard input.\n"
  35 "PATTERN, by default, is a perl compatible regular expression.\n\n"
  36 "Options:\n"
  37 "  -h, --help              print this message\n"
  38 "  -v, --invert-match      select nonmatching entries\n"
  39 "  -c, --count             output only the number of matching sequences\n"
  40     );
  41 }
  42
  43 static int invert_flag;
  44 static int help_flag;
  45 static int count_flag;
  46
  47
  48
  49 void print_fastq_entry( FILE* fout, kseq_t* seq )
  50 {
  51     fprintf(fout, "@%s\n%s\n+%s\n%s\n",
  52                   seq->name.s,
  53                   seq->seq.s,
  54                   seq->comment.s,
  55                   seq->qual.s );
  56 }
  57
  58
  59 void fastq_grep( gzFile fin, FILE* fout, pcre* re )
  60 {
  61     int rc;
  62     int ovector[3];
  63     size_t count = 0;
  64
  65     kseq_t* seq;
  66     seq = kseq_init(fin);
  67
  68     while (kseq_read(seq) >= 0) {
  69         rc = pcre_exec(re,          /* pattern */
  70                        NULL,        /* extre data */
  71                        seq->seq.s,  /* subject */
  72                        seq->seq.l,  /* subject length */
  73                        0,           /* subject offset */
  74                        0,           /* options */
  75                        ovector,     /* output vector */
  76                        3         ); /* output vector length */
  77
  78         if ((invert_flag && rc == PCRE_ERROR_NOMATCH) || rc >= 0) {
  79             if (count_flag) count++;
  80             else            print_fastq_entry(fout, seq);
  81         }
  82     }
  83
  84     kseq_destroy(seq);
  85
  86     if (count_flag) fprintf(fout, "%zu\n", count);
  87 }
  88
  89
  90
  91 int main(int argc, char* argv[])
  92 {
  93     SET_BINARY_MODE(stdin);
  94     SET_BINARY_MODE(stdout);
  95
  96     const char* pat;
  97     pcre* re;
  98     const char* pat_error;
  99     int pat_error_offset;
 100
 101     FILE*  fin;
 102     gzFile gzfin;
 103
 104
 105     invert_flag  = 0;
 106     help_flag    = 0;
 107     count_flag   = 0;
 108
 109     int opt;
 110     int opt_idx;
 111
 112
 113     static struct option long_options[] =
 114         {
 115           {"help", no_argument, &help_flag, 1},
 116           {"invert-match", no_argument, &invert_flag, 1},
 117           {"count", no_argument, &count_flag, 1},
 118           {0, 0, 0, 0}
 119         };
 120
 121     while (1) {
 122         opt = getopt_long(argc, argv, "hvc", long_options, &opt_idx);
 123
 124         if( opt == -1 ) break;
 125
 126         switch (opt) {
 127             case 0:
 128                 if (long_options[opt_idx].flag != 0) break;
 129                 if (optarg) {
 130                 }
 131                 break;
 132
 133             case 'h':
 134                 help_flag = 1;
 135                 break;
 136
 137             case 'v':
 138                 invert_flag = 1;
 139                 break;
 140
 141             case 'c':
 142                 count_flag = 1;
 143                 break;
 144
 145             case '?':
 146                 return 1;
 147
 148             default:
 149                 abort();
 150         }
 151     }
 152
 153     if (help_flag) {
 154         print_help();
 155         return 0;
 156     }
 157
 158     if (optind >= argc) {
 159         fprintf(stderr, "A pattern must be specified.\n");
 160         return 1;
 161     }
 162
 163     pat = argv[optind++];
 164     re = pcre_compile( pat, PCRE_CASELESS, &pat_error, &pat_error_offset, NULL );
 165
 166
 167     if (re == NULL) {
 168         fprintf(stderr, "Syntax error in PCRE pattern at offset: %d: %s\n",
 169                 pat_error_offset, pat_error );
 170         return 1;
 171     }
 172
 173
 174     if (optind >= argc || (argc - optind == 1 && strcmp(argv[optind],"-") == 0)) {
 175         gzfin = gzdopen( fileno(stdin), "rb" );
 176         if (gzfin == NULL) {
 177             fprintf(stderr, "Malformed file 'stdin'.\n");
 178             return 1;
 179         }
 180
 181         fastq_grep(gzfin, stdout, re);
 182
 183         gzclose(gzfin);
 184     }
 185     else {
 186         for (; optind < argc; optind++) {
 187             fin = fopen(argv[optind], "rb");
 188             if (fin == NULL) {
 189                 fprintf(stderr, "No such file '%s'.\n", argv[optind]);
 190                 continue;
 191             }
 192
 193             gzfin = gzdopen(fileno(fin), "rb");
 194             if (gzfin == NULL) {
 195                 fprintf(stderr, "Malformed file '%s'.\n", argv[optind]);
 196                 continue;
 197             }
 198
 199             fastq_grep(gzfin, stdout, re);
 200
 201             gzclose(gzfin);
 202         }
 203     }
 204
 205     pcre_free(re);
 206
 207     return 0;
 208 }
 209
 210
 211