main.c

   1 #include <string.h>
   2 #include <unistd.h>
   3 #include <stdlib.h>
   4 #include <stdio.h>
   5 #include <sys/stat.h>
   6 #include <errno.h>
   7 #include "bgzf.h"
   8 #include "tabix.h"
   9 #include "knetfile.h"
  10
  11 #define PACKAGE_VERSION "0.2.5 (r1005)"
  12
  13 #define error(...) { fprintf(stderr,__VA_ARGS__); return -1; }
  14
  15 int reheader_file(const char *header, const char *file, int meta)
  16 {
  17     BGZF *fp = bgzf_open(file,"r");
  18     if (bgzf_read_block(fp) != 0 || !fp->block_length)
  19         return -1;
  20
  21     char *buffer = fp->uncompressed_block;
  22     int skip_until = 0;
  23
  24     if ( buffer[0]==meta )
  25     {
  26         skip_until = 1;
  27
  28         // Skip the header
  29         while (1)
  30         {
  31             if ( buffer[skip_until]=='\n' )
  32             {
  33                 skip_until++;
  34                 if ( skip_until>=fp->block_length )
  35                 {
  36                     if (bgzf_read_block(fp) != 0 || !fp->block_length)
  37                         error("no body?\n");
  38                     skip_until = 0;
  39                 }
  40                 // The header has finished
  41                 if ( buffer[skip_until]!=meta ) break;
  42             }
  43             skip_until++;
  44             if ( skip_until>=fp->block_length )
  45             {
  46                 if (bgzf_read_block(fp) != 0 || !fp->block_length)
  47                     error("no body?\n");
  48                 skip_until = 0;
  49             }
  50         }
  51     }
  52
  53     FILE *fh = fopen(header,"r");
  54     if ( !fh )
  55         error("%s: %s", header,strerror(errno));
  56     int page_size = getpagesize();
  57     char *buf = valloc(page_size);
  58     BGZF *bgzf_out = bgzf_dopen(fileno(stdout), "w");
  59     ssize_t nread;
  60     while ( (nread=fread(buf,1,page_size-1,fh))>0 )
  61     {
  62         if ( nread<page_size-1 && buf[nread-1]!='\n' )
  63             buf[nread++] = '\n';
  64         if (bgzf_write(bgzf_out, buf, nread) < 0) error("Error: %d\n",bgzf_out->errcode);
  65     }
  66     fclose(fh);
  67
  68     if ( fp->block_length - skip_until > 0 )
  69     {
  70         if (bgzf_write(bgzf_out, buffer+skip_until, fp->block_length-skip_until) < 0)
  71             error("Error: %d\n",fp->errcode);
  72     }
  73     if (bgzf_flush(bgzf_out) < 0)
  74         error("Error: %d\n",bgzf_out->errcode);
  75
  76     while (1)
  77     {
  78 #ifdef _USE_KNETFILE
  79         nread = knet_read(fp->fp, buf, page_size);
  80 #else
  81         nread = fread(buf, 1, page_size, fp->fp);
  82 #endif
  83         if ( nread<=0 )
  84             break;
  85
  86         int count = fwrite(buf, 1, nread, bgzf_out->fp);
  87         if (count != nread)
  88             error("Write failed, wrote %d instead of %d bytes.\n", count,(int)nread);
  89     }
  90
  91     if (bgzf_close(bgzf_out) < 0)
  92         error("Error: %d\n",bgzf_out->errcode);
  93
  94     return 0;
  95 }
  96
  97
  98 int main(int argc, char *argv[])
  99 {
 100         int c, skip = -1, meta = -1, list_chrms = 0, force = 0, print_header = 0, print_only_header = 0, bed_reg = 0;
 101         ti_conf_t conf = ti_conf_gff, *conf_ptr = NULL;
 102     const char *reheader = NULL;
 103         while ((c = getopt(argc, argv, "p:s:b:e:0S:c:lhHfBr:")) >= 0) {
 104                 switch (c) {
 105                 case 'B': bed_reg = 1; break;
 106                 case '0': conf.preset |= TI_FLAG_UCSC; break;
 107                 case 'S': skip = atoi(optarg); break;
 108                 case 'c': meta = optarg[0]; break;
 109                 case 'p':
 110                         if (strcmp(optarg, "gff") == 0) conf_ptr = &ti_conf_gff;
 111                         else if (strcmp(optarg, "bed") == 0) conf_ptr = &ti_conf_bed;
 112                         else if (strcmp(optarg, "sam") == 0) conf_ptr = &ti_conf_sam;
 113                         else if (strcmp(optarg, "vcf") == 0 || strcmp(optarg, "vcf4") == 0) conf_ptr = &ti_conf_vcf;
 114                         else if (strcmp(optarg, "psltbl") == 0) conf_ptr = &ti_conf_psltbl;
 115                         else {
 116                                 fprintf(stderr, "[main] unrecognized preset '%s'\n", optarg);
 117                                 return 1;
 118                         }
 119                         break;
 120                 case 's': conf.sc = atoi(optarg); break;
 121                 case 'b': conf.bc = atoi(optarg); break;
 122                 case 'e': conf.ec = atoi(optarg); break;
 123         case 'l': list_chrms = 1; break;
 124         case 'h': print_header = 1; break;
 125         case 'H': print_only_header = 1; break;
 126                 case 'f': force = 1; break;
 127         case 'r': reheader = optarg; break;
 128                 }
 129         }
 130         if (optind == argc) {
 131                 fprintf(stderr, "\n");
 132                 fprintf(stderr, "Program: tabix (TAB-delimited file InderXer)\n");
 133                 fprintf(stderr, "Version: %s\n\n", PACKAGE_VERSION);
 134                 fprintf(stderr, "Usage:   tabix <in.tab.bgz> [region1 [region2 [...]]]\n\n");
 135                 fprintf(stderr, "Options: -p STR     preset: gff, bed, sam, vcf, psltbl [gff]\n");
 136                 fprintf(stderr, "         -s INT     sequence name column [1]\n");
 137                 fprintf(stderr, "         -b INT     start column [4]\n");
 138                 fprintf(stderr, "         -e INT     end column; can be identical to '-b' [5]\n");
 139                 fprintf(stderr, "         -S INT     skip first INT lines [0]\n");
 140                 fprintf(stderr, "         -c CHAR    symbol for comment/meta lines [#]\n");
 141             fprintf(stderr, "         -r FILE    replace the header with the content of FILE [null]\n");
 142                 fprintf(stderr, "         -B         region1 is a BED file (entire file will be read)\n");
 143                 fprintf(stderr, "         -0         zero-based coordinate\n");
 144                 fprintf(stderr, "         -h         print also the header lines\n");
 145                 fprintf(stderr, "         -H         print only the header lines\n");
 146                 fprintf(stderr, "         -l         list chromosome names\n");
 147                 fprintf(stderr, "         -f         force to overwrite the index\n");
 148                 fprintf(stderr, "\n");
 149                 return 1;
 150         }
 151     if ( !conf_ptr )
 152     {
 153         int l = strlen(argv[optind]);
 154         int strcasecmp(const char *s1, const char *s2);
 155         if (l>=7 && strcasecmp(argv[optind]+l-7, ".gff.gz") == 0) conf_ptr = &ti_conf_gff;
 156         else if (l>=7 && strcasecmp(argv[optind]+l-7, ".bed.gz") == 0) conf_ptr = &ti_conf_bed;
 157         else if (l>=7 && strcasecmp(argv[optind]+l-7, ".sam.gz") == 0) conf_ptr = &ti_conf_sam;
 158         else if (l>=7 && strcasecmp(argv[optind]+l-7, ".vcf.gz") == 0) conf_ptr = &ti_conf_vcf;
 159         else if (l>=10 && strcasecmp(argv[optind]+l-10, ".psltbl.gz") == 0) conf_ptr = &ti_conf_psltbl;
 160     }
 161     if ( conf_ptr )
 162         conf = *conf_ptr;
 163
 164         if (skip >= 0) conf.line_skip = skip;
 165         if (meta >= 0) conf.meta_char = meta;
 166     if (list_chrms) {
 167                 ti_index_t *idx;
 168                 int i, n;
 169                 const char **names;
 170                 idx = ti_index_load(argv[optind]);
 171                 if (idx == 0) {
 172                         fprintf(stderr, "[main] fail to load the index file.\n");
 173                         return 1;
 174                 }
 175                 names = ti_seqname(idx, &n);
 176                 for (i = 0; i < n; ++i) printf("%s\n", names[i]);
 177                 free(names);
 178                 ti_index_destroy(idx);
 179                 return 0;
 180         }
 181     if (reheader)
 182         return reheader_file(reheader,argv[optind],conf.meta_char);
 183
 184         struct stat stat_tbi,stat_vcf;
 185     char *fnidx = calloc(strlen(argv[optind]) + 5, 1);
 186         strcat(strcpy(fnidx, argv[optind]), ".tbi");
 187
 188         if (optind + 1 == argc && !print_only_header) {
 189                 if (force == 0) {
 190                         if (stat(fnidx, &stat_tbi) == 0)
 191             {
 192                 // Before complaining, check if the VCF file isn't newer. This is a common source of errors,
 193                 //  people tend not to notice that tabix failed
 194                 stat(argv[optind], &stat_vcf);
 195                 if ( stat_vcf.st_mtime <= stat_tbi.st_mtime )
 196                 {
 197                     fprintf(stderr, "[tabix] the index file exists. Please use '-f' to overwrite.\n");
 198                     free(fnidx);
 199                     return 1;
 200                 }
 201                         }
 202                 }
 203         if ( bgzf_is_bgzf(argv[optind])!=1 )
 204         {
 205             fprintf(stderr,"[tabix] was bgzip used to compress this file? %s\n", argv[optind]);
 206             free(fnidx);
 207             return 1;
 208         }
 209         if ( !conf_ptr )
 210         {
 211             // Building the index but the file type was neither recognised nor given. If no custom change
 212             //  has been made, warn the user that GFF is used
 213             if ( conf.preset==ti_conf_gff.preset
 214                 && conf.sc==ti_conf_gff.sc
 215                 && conf.bc==ti_conf_gff.bc
 216                 && conf.ec==ti_conf_gff.ec
 217                 && conf.meta_char==ti_conf_gff.meta_char
 218                 && conf.line_skip==ti_conf_gff.line_skip )
 219                 fprintf(stderr,"[tabix] The file type not recognised and -p not given, using the preset [gff].\n");
 220         }
 221                 return ti_index_build(argv[optind], &conf);
 222         }
 223         { // retrieve
 224                 tabix_t *t;
 225         // On some systems, stat on non-existent files returns undefined value for sm_mtime, the user had to use -f
 226         int is_remote = (strstr(fnidx, "ftp://") == fnidx || strstr(fnidx, "http://") == fnidx) ? 1 : 0;
 227         if ( !is_remote )
 228         {
 229             // Common source of errors: new VCF is used with an old index
 230             stat(fnidx, &stat_tbi);
 231             stat(argv[optind], &stat_vcf);
 232             if ( force==0 && stat_vcf.st_mtime > stat_tbi.st_mtime )
 233             {
 234                 fprintf(stderr, "[tabix] the index file either does not exist or is older than the vcf file. Please reindex.\n");
 235                 free(fnidx);
 236                 return 1;
 237             }
 238         }
 239         free(fnidx);
 240
 241                 if ((t = ti_open(argv[optind], 0)) == 0) {
 242                         fprintf(stderr, "[main] fail to open the data file.\n");
 243                         return 1;
 244                 }
 245         if ( print_only_header )
 246         {
 247             ti_iter_t iter;
 248             const char *s;
 249             int len;
 250             if (ti_lazy_index_load(t) < 0 && bed_reg == 0) {
 251                 fprintf(stderr,"[tabix] failed to load the index file.\n");
 252                 return 1;
 253             }
 254             const ti_conf_t *idxconf = ti_get_conf(t->idx);
 255             iter = ti_query(t, 0, 0, 0);
 256             while ((s = ti_read(t, iter, &len)) != 0) {
 257                 if ((int)(*s) != idxconf->meta_char) break;
 258                 fputs(s, stdout); fputc('\n', stdout);
 259             }
 260             ti_iter_destroy(iter);
 261             return 0;
 262         }
 263
 264                 if (strcmp(argv[optind+1], ".") == 0) { // retrieve all
 265                         ti_iter_t iter;
 266                         const char *s;
 267                         int len;
 268                         iter = ti_query(t, 0, 0, 0);
 269                         while ((s = ti_read(t, iter, &len)) != 0) {
 270                                 fputs(s, stdout); fputc('\n', stdout);
 271                         }
 272                         ti_iter_destroy(iter);
 273                 } else { // retrieve from specified regions
 274                         int i, len;
 275             ti_iter_t iter;
 276             const char *s;
 277                         const ti_conf_t *idxconf;
 278
 279                         if (ti_lazy_index_load(t) < 0 && bed_reg == 0) {
 280                 fprintf(stderr,"[tabix] failed to load the index file.\n");
 281                 return 1;
 282             }
 283                         idxconf = ti_get_conf(t->idx);
 284
 285             if ( print_header )
 286             {
 287                 // If requested, print the header lines here
 288                 iter = ti_query(t, 0, 0, 0);
 289                 while ((s = ti_read(t, iter, &len)) != 0) {
 290                     if ((int)(*s) != idxconf->meta_char) break;
 291                     fputs(s, stdout); fputc('\n', stdout);
 292                 }
 293                 ti_iter_destroy(iter);
 294             }
 295                         if (bed_reg) {
 296                                 extern int bed_overlap(const void *_h, const char *chr, int beg, int end);
 297                                 extern void *bed_read(const char *fn);
 298                                 extern void bed_destroy(void *_h);
 299
 300                                 const ti_conf_t *conf_ = idxconf? idxconf : &conf; // use the index file if available
 301                                 void *bed = bed_read(argv[optind+1]); // load the BED file
 302                                 ti_interval_t intv;
 303
 304                                 if (bed == 0) {
 305                                         fprintf(stderr, "[main] fail to read the BED file.\n");
 306                                         return 1;
 307                                 }
 308                                 iter = ti_query(t, 0, 0, 0);
 309                                 while ((s = ti_read(t, iter, &len)) != 0) {
 310                                         int c;
 311                                         ti_get_intv(conf_, len, (char*)s, &intv);
 312                                         c = *intv.se; *intv.se = '\0';
 313                                         if (bed_overlap(bed, intv.ss, intv.beg, intv.end)) {
 314                                                 *intv.se = c;
 315                                                 puts(s);
 316                                         }
 317                                         *intv.se = c;
 318                                 }
 319                 ti_iter_destroy(iter);
 320                                 bed_destroy(bed);
 321                         } else {
 322                                 for (i = optind + 1; i < argc; ++i) {
 323                                         int tid, beg, end;
 324                                         if (ti_parse_region(t->idx, argv[i], &tid, &beg, &end) == 0) {
 325                                                 iter = ti_queryi(t, tid, beg, end);
 326                                                         while ((s = ti_read(t, iter, &len)) != 0) {
 327                                                         fputs(s, stdout); fputc('\n', stdout);
 328                                                 }
 329                                                 ti_iter_destroy(iter);
 330                                         }
 331                     // else fprintf(stderr, "[main] invalid region: unknown target name or minus interval.\n");
 332                                 }
 333                         }
 334                 }
 335                 ti_close(t);
 336         }
 337         return 0;
 338 }