src/fastq-match.c

   1
   2 /* Smith-Waterman alignments against sequences within a fastq file.
   3  *
   4  * Daniel Jones <dcjones@cs.washington.edu>
   5  */
   6
   7 #include <stdlib.h>
   8 #include <zlib.h>
   9 #include <getopt.h>
  10
  11 #include "swsse2/blosum62.h"
  12 #include "swsse2/swsse2.h"
  13 #include "swsse2/matrix.h"
  14 #include "swsse2/swstriped.h"
  15 #include "kseq.h"
  16 KSEQ_INIT(gzFile, gzread)
  17
  18
  19 #if defined(MSDOS) || defined(OS2) || defined(WIN32) || defined(__CYGWIN__)
  20 #  include <fcntl.h>
  21 #  include <io.h>
  22 #  define SET_BINARY_MODE(file) setmode(fileno(file), O_BINARY)
  23 #else
  24 #  define SET_BINARY_MODE(file)
  25 #endif
  26
  27
  28 static int help_flag;
  29
  30
  31 void print_help()
  32 {
  33     fprintf( stderr,
  34 "fastq-match [OPTION]... QUERY [FILE]...\n"
  35 "Perform Smith-Waterman local alignment of a query sequence\n"
  36 "against each sequence in a fastq file.\n\n"
  37 "Options:\n"
  38 "  -h, --help              print this message\n"
  39     );
  40 }
  41
  42
  43 void convert_sequence(unsigned char* s, int n)
  44 {
  45     int i;
  46     for (i = 0; i < n; i++) {
  47         s[i] = (char)AMINO_ACID_VALUE[(int)s[i]];
  48     }
  49 }
  50
  51
  52 void fastq_match(gzFile fin, FILE* fout,
  53                  SwStripedData* sw_data,
  54                  unsigned char* query, int n,
  55                  SEARCH_OPTIONS* options)
  56 {
  57     int score;
  58
  59     int gap_init  = -(options->gapInit + options->gapExt);
  60     int gap_ext   = -options->gapExt;
  61     int threshold = options->threshold;
  62
  63     kseq_t* seq;
  64     seq = kseq_init(fin);
  65
  66     while (kseq_read(seq) >= 0) {
  67         fprintf(fout, "%s\t", seq->seq.s);
  68
  69         convert_sequence((unsigned char*)seq->seq.s, seq->seq.l);
  70
  71         score = swStripedByte(query, n,
  72                               (unsigned char*)seq->seq.s, seq->seq.l,
  73                               gap_init, gap_ext,
  74                               sw_data->pvbQueryProf,
  75                               sw_data->pvH1,
  76                               sw_data->pvH2,
  77                               sw_data->pvE,
  78                               sw_data->bias);
  79         if (score >= 255) {
  80             score = swStripedWord(query, n,
  81                                   (unsigned char*)seq->seq.s, seq->seq.l,
  82                                   gap_init, gap_ext,
  83                                   sw_data->pvbQueryProf,
  84                                   sw_data->pvH1,
  85                                   sw_data->pvH2,
  86                                   sw_data->pvE);
  87         }
  88
  89         fprintf(fout, "%d\n", score);
  90     }
  91
  92     kseq_destroy(seq);
  93 }
  94
  95
  96
  97 int main(int argc, char* argv[])
  98 {
  99     SET_BINARY_MODE(stdin);
 100     SET_BINARY_MODE(stdout);
 101
 102     unsigned char* query;
 103     int query_len;
 104     SwStripedData* sw_data;
 105     signed char* mat = blosum62;
 106     SEARCH_OPTIONS options;
 107
 108     options.gapInit   = -10;
 109     options.gapExt    = -2;
 110     options.threshold = -1;
 111
 112     FILE*  fin;
 113     gzFile gzfin;
 114
 115     help_flag = 0;
 116
 117     int opt;
 118     int opt_idx;
 119
 120     static struct option long_options[] =
 121         {
 122           {"help", no_argument, &help_flag, 1},
 123           {0, 0, 0, 0}
 124         };
 125
 126
 127     while (1) {
 128         opt = getopt_long(argc, argv, "h", long_options, &opt_idx);
 129
 130         if( opt == -1 ) break;
 131
 132         switch (opt) {
 133             case 0:
 134                 if (long_options[opt_idx].flag != 0) break;
 135                 if (optarg) {
 136                 }
 137                 break;
 138
 139             case 'h':
 140                 help_flag = 1;
 141                 break;
 142
 143             case '?':
 144                 return 1;
 145
 146             default:
 147                 abort();
 148         }
 149     }
 150
 151     if (help_flag) {
 152         print_help();
 153         return 0;
 154     }
 155
 156     if (optind >= argc) {
 157         fprintf(stderr, "A query sequence must be specified.\n");
 158         return 1;
 159     }
 160
 161     query = (unsigned char*)argv[optind++];
 162     query_len = strlen((char*)query);
 163     convert_sequence(query, query_len);
 164
 165     sw_data = swStripedInit(query, query_len, mat);
 166
 167     if (optind >= argc || (argc - optind == 1 && strcmp(argv[optind],"-") == 0)) {
 168         gzfin = gzdopen( fileno(stdin), "rb" );
 169         if (gzfin == NULL) {
 170             fprintf(stderr, "Malformed file 'stdin'.\n");
 171             return 1;
 172         }
 173
 174         fastq_match(gzfin, stdout, sw_data, query, query_len, &options);
 175
 176         gzclose(gzfin);
 177     }
 178     else {
 179         for (; optind < argc; optind++) {
 180             fin = fopen(argv[optind], "rb");
 181             if (fin == NULL) {
 182                 fprintf(stderr, "No such file '%s'.\n", argv[optind]);
 183                 continue;
 184             }
 185
 186             gzfin = gzdopen(fileno(fin), "rb");
 187             if (gzfin == NULL) {
 188                 fprintf(stderr, "Malformed file '%s'.\n", argv[optind]);
 189                 continue;
 190             }
 191
 192             fastq_match(gzfin, stdout, sw_data, query, query_len, &options);
 193
 194             gzclose(gzfin);
 195         }
 196     }
 197
 198     swStripedComplete(sw_data);
 199
 200     return 0;
 201 }
 202
 203
 204