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     );
  40 }
  41
  42 static int invert_flag;
  43 static int help_flag;
  44 static int zout_flag;
  45
  46
  47
  48 void print_fastq_entry( FILE* fout, kseq_t* seq )
  49 {
  50     fprintf(fout, "@%s\n%s\n+%s\n%s\n",
  51                   seq->name.s,
  52                   seq->seq.s,
  53                   seq->comment.s,
  54                   seq->qual.s );
  55 }
  56
  57
  58 void fastq_grep( gzFile fin, FILE* fout, pcre* re )
  59 {
  60     int rc;
  61     int ovector[3];
  62
  63     kseq_t* seq;
  64     seq = kseq_init(fin);
  65
  66     while (kseq_read(seq) >= 0) {
  67         rc = pcre_exec(re,          /* pattern */
  68                        NULL,        /* extre data */
  69                        seq->seq.s,  /* subject */
  70                        seq->seq.l,  /* subject length */
  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) || rc >= 0) {
  77             print_fastq_entry( fout, seq );
  78         }
  79     }
  80
  81     kseq_destroy(seq);
  82 }
  83
  84
  85
  86 int main(int argc, char* argv[])
  87 {
  88     const char* pat;
  89     pcre* re;
  90     const char* pat_error;
  91     int pat_error_offset;
  92
  93     FILE*  fin;
  94     gzFile gzfin;
  95
  96
  97     invert_flag = 0;
  98     help_flag   = 0;
  99     zout_flag   = 0;
 100
 101     int opt;
 102     int opt_idx;
 103
 104
 105     static struct option long_options[] =
 106         {
 107           {"help", no_argument, &help_flag, 1},
 108           {"invert-match", no_argument, &invert_flag, 1},
 109           {0, 0, 0, 0}
 110         };
 111
 112     while (1) {
 113         opt = getopt_long(argc, argv, "hv", long_options, &opt_idx);
 114
 115         if( opt == -1 ) break;
 116
 117         switch (opt) {
 118             case 0:
 119                 if (long_options[opt_idx].flag != 0) break;
 120                 if (optarg) {
 121                     /* TODO */
 122                 }
 123                 break;
 124
 125             case 'h':
 126                 help_flag = 1;
 127                 break;
 128
 129             case 'v':
 130                 invert_flag = 1;
 131                 break;
 132
 133             default:
 134                 abort();
 135         }
 136     }
 137
 138     if (help_flag) {
 139         print_help();
 140         return 0;
 141     }
 142
 143     if (optind >= argc) {
 144         fprintf(stderr, "A pattern must be specified.\n");
 145         return 1;
 146     }
 147
 148     pat = argv[optind++];
 149     re = pcre_compile( pat, PCRE_CASELESS, &pat_error, &pat_error_offset, NULL );
 150
 151
 152     if (re == NULL) {
 153         fprintf(stderr, "Syntax error in PCRE pattern at offset: %d: %s\n",
 154                 pat_error_offset, pat_error );
 155         return 1;
 156     }
 157
 158
 159     if (optind >= argc || (argc - optind == 1 && strcmp(argv[optind],"-") == 0)) {
 160         gzfin = gzdopen( fileno(stdin), "rb" );
 161         if (gzfin == NULL) {
 162             fprintf(stderr, "Malformed file 'stdin'.\n");
 163             return 1;
 164         }
 165
 166         fastq_grep(gzfin, stdout, re);
 167
 168         gzclose(gzfin);
 169     }
 170     else {
 171         for (; optind < argc; optind++) {
 172             fin = fopen(argv[optind], "rb");
 173             if (fin == NULL) {
 174                 fprintf(stderr, "No such file '%s'.\n", argv[optind]);
 175                 continue;
 176             }
 177
 178             gzfin = gzdopen(fileno(fin), "rb");
 179             if (gzfin == NULL) {
 180                 fprintf(stderr, "Malformed file '%s'.\n", argv[optind]);
 181                 continue;
 182             }
 183
 184             fastq_grep(gzfin, stdout, re);
 185
 186             gzclose(gzfin);
 187         }
 188     }
 189
 190     pcre_free(re);
 191
 192     return 0;
 193 }
 194
 195
 196