faidx.c

   1 #include <ctype.h>
   2 #include <string.h>
   3 #include <stdlib.h>
   4 #include <stdio.h>
   5 #include <stdint.h>
   6 #include "faidx.h"
   7 #include "khash.h"
   8
   9 typedef struct {
  10         uint64_t len:32, line_len:16, line_blen:16;
  11         uint64_t offset;
  12 } faidx1_t;
  13 KHASH_MAP_INIT_STR(s, faidx1_t)
  14
  15 #ifndef _NO_RAZF
  16 #include "razf.h"
  17 #else
  18 #ifdef _WIN32
  19 #define ftello(fp) ftell(fp)
  20 #define fseeko(fp, offset, whence) fseek(fp, offset, whence)
  21 #else
  22 extern off_t ftello(FILE *stream);
  23 extern int fseeko(FILE *stream, off_t offset, int whence);
  24 #endif
  25 #define RAZF FILE
  26 #define razf_read(fp, buf, size) fread(buf, 1, size, fp)
  27 #define razf_open(fn, mode) fopen(fn, mode)
  28 #define razf_close(fp) fclose(fp)
  29 #define razf_seek(fp, offset, whence) fseeko(fp, offset, whence)
  30 #define razf_tell(fp) ftello(fp)
  31 #endif
  32 #ifdef _USE_KNETFILE
  33 #include "knetfile.h"
  34 #endif
  35
  36 struct __faidx_t {
  37         RAZF *rz;
  38         int n, m;
  39         char **name;
  40         khash_t(s) *hash;
  41 };
  42
  43 #ifndef kroundup32
  44 #define kroundup32(x) (--(x), (x)|=(x)>>1, (x)|=(x)>>2, (x)|=(x)>>4, (x)|=(x)>>8, (x)|=(x)>>16, ++(x))
  45 #endif
  46
  47 static inline void fai_insert_index(faidx_t *idx, const char *name, int len, int line_len, int line_blen, uint64_t offset)
  48 {
  49         khint_t k;
  50         int ret;
  51         faidx1_t t;
  52         if (idx->n == idx->m) {
  53                 idx->m = idx->m? idx->m<<1 : 16;
  54                 idx->name = (char**)realloc(idx->name, sizeof(void*) * idx->m);
  55         }
  56         idx->name[idx->n] = strdup(name);
  57         k = kh_put(s, idx->hash, idx->name[idx->n], &ret);
  58         t.len = len; t.line_len = line_len; t.line_blen = line_blen; t.offset = offset;
  59         kh_value(idx->hash, k) = t;
  60         ++idx->n;
  61 }
  62
  63 faidx_t *fai_build_core(RAZF *rz)
  64 {
  65         char c, *name;
  66         int l_name, m_name, ret;
  67         int len, line_len, line_blen, state;
  68         int l1, l2;
  69         faidx_t *idx;
  70         uint64_t offset;
  71
  72         idx = (faidx_t*)calloc(1, sizeof(faidx_t));
  73         idx->hash = kh_init(s);
  74         name = 0; l_name = m_name = 0;
  75         len = line_len = line_blen = -1; state = 0; l1 = l2 = -1; offset = 0;
  76         while (razf_read(rz, &c, 1)) {
  77                 if (c == '\n') { // an empty line
  78                         if (state == 1) {
  79                                 offset = razf_tell(rz);
  80                                 continue;
  81                         } else if ((state == 0 && len < 0) || state == 2) continue;
  82                 }
  83                 if (c == '>') { // fasta header
  84                         if (len >= 0)
  85                                 fai_insert_index(idx, name, len, line_len, line_blen, offset);
  86                         l_name = 0;
  87                         while ((ret = razf_read(rz, &c, 1)) != 0 && !isspace(c)) {
  88                                 if (m_name < l_name + 2) {
  89                                         m_name = l_name + 2;
  90                                         kroundup32(m_name);
  91                                         name = (char*)realloc(name, m_name);
  92                                 }
  93                                 name[l_name++] = c;
  94                         }
  95                         name[l_name] = '\0';
  96                         if (ret == 0) {
  97                                 fprintf(stderr, "[fai_build_core] the last entry has no sequence\n");
  98                                 free(name); fai_destroy(idx);
  99                                 return 0;
 100                         }
 101                         if (c != '\n') while (razf_read(rz, &c, 1) && c != '\n');
 102                         state = 1; len = 0;
 103                         offset = razf_tell(rz);
 104                 } else {
 105                         if (state == 3) {
 106                                 fprintf(stderr, "[fai_build_core] inlined empty line is not allowed in sequence '%s'.\n", name);
 107                                 free(name); fai_destroy(idx);
 108                                 return 0;
 109                         }
 110                         if (state == 2) state = 3;
 111                         l1 = l2 = 0;
 112                         do {
 113                                 ++l1;
 114                                 if (isgraph(c)) ++l2;
 115                         } while ((ret = razf_read(rz, &c, 1)) && c != '\n');
 116                         if (state == 3 && l2) {
 117                                 fprintf(stderr, "[fai_build_core] different line length in sequence '%s'.\n", name);
 118                                 free(name); fai_destroy(idx);
 119                                 return 0;
 120                         }
 121                         ++l1; len += l2;
 122                         if (l2 >= 0x10000) {
 123                                 fprintf(stderr, "[fai_build_core] line length exceeds 65535 in sequence '%s'.\n", name);
 124                                 free(name); fai_destroy(idx);
 125                                 return 0;
 126                         }
 127                         if (state == 1) line_len = l1, line_blen = l2, state = 0;
 128                         else if (state == 0) {
 129                                 if (l1 != line_len || l2 != line_blen) state = 2;
 130                         }
 131                 }
 132         }
 133         fai_insert_index(idx, name, len, line_len, line_blen, offset);
 134         free(name);
 135         return idx;
 136 }
 137
 138 void fai_save(const faidx_t *fai, FILE *fp)
 139 {
 140         khint_t k;
 141         int i;
 142         for (i = 0; i < fai->n; ++i) {
 143                 faidx1_t x;
 144                 k = kh_get(s, fai->hash, fai->name[i]);
 145                 x = kh_value(fai->hash, k);
 146 #ifdef _WIN32
 147                 fprintf(fp, "%s\t%d\t%ld\t%d\t%d\n", fai->name[i], (int)x.len, (long)x.offset, (int)x.line_blen, (int)x.line_len);
 148 #else
 149                 fprintf(fp, "%s\t%d\t%lld\t%d\t%d\n", fai->name[i], (int)x.len, (long long)x.offset, (int)x.line_blen, (int)x.line_len);
 150 #endif
 151         }
 152 }
 153
 154 faidx_t *fai_read(FILE *fp)
 155 {
 156         faidx_t *fai;
 157         char *buf, *p;
 158         int len, line_len, line_blen;
 159 #ifdef _WIN32
 160         long offset;
 161 #else
 162         long long offset;
 163 #endif
 164         fai = (faidx_t*)calloc(1, sizeof(faidx_t));
 165         fai->hash = kh_init(s);
 166         buf = (char*)calloc(0x10000, 1);
 167         while (!feof(fp) && fgets(buf, 0x10000, fp)) {
 168                 for (p = buf; *p && isgraph(*p); ++p);
 169                 *p = 0; ++p;
 170 #ifdef _WIN32
 171                 sscanf(p, "%d%ld%d%d", &len, &offset, &line_blen, &line_len);
 172 #else
 173                 sscanf(p, "%d%lld%d%d", &len, &offset, &line_blen, &line_len);
 174 #endif
 175                 fai_insert_index(fai, buf, len, line_len, line_blen, offset);
 176         }
 177         free(buf);
 178         return fai;
 179 }
 180
 181 void fai_destroy(faidx_t *fai)
 182 {
 183         int i;
 184         for (i = 0; i < fai->n; ++i) free(fai->name[i]);
 185         free(fai->name);
 186         kh_destroy(s, fai->hash);
 187         if (fai->rz) razf_close(fai->rz);
 188         free(fai);
 189 }
 190
 191 int fai_build(const char *fn)
 192 {
 193         char *str;
 194         RAZF *rz;
 195         FILE *fp;
 196         faidx_t *fai;
 197         str = (char*)calloc(strlen(fn) + 5, 1);
 198         sprintf(str, "%s.fai", fn);
 199         rz = razf_open(fn, "r");
 200         if (rz == 0) {
 201                 fprintf(stderr, "[fai_build] fail to open the FASTA file %s\n",fn);
 202                 free(str);
 203                 return -1;
 204         }
 205         fai = fai_build_core(rz);
 206         razf_close(rz);
 207         fp = fopen(str, "wb");
 208         if (fp == 0) {
 209                 fprintf(stderr, "[fai_build] fail to write FASTA index %s\n",str);
 210                 fai_destroy(fai); free(str);
 211                 return -1;
 212         }
 213         fai_save(fai, fp);
 214         fclose(fp);
 215         free(str);
 216         fai_destroy(fai);
 217         return 0;
 218 }
 219
 220 #ifdef _USE_KNETFILE
 221 FILE *download_and_open(const char *fn)
 222 {
 223     const int buf_size = 1 * 1024 * 1024;
 224     uint8_t *buf;
 225     FILE *fp;
 226     knetFile *fp_remote;
 227     const char *url = fn;
 228     const char *p;
 229     int l = strlen(fn);
 230     for (p = fn + l - 1; p >= fn; --p)
 231         if (*p == '/') break;
 232     fn = p + 1;
 233
 234     // First try to open a local copy
 235     fp = fopen(fn, "r");
 236     if (fp)
 237         return fp;
 238
 239     // If failed, download from remote and open
 240     fp_remote = knet_open(url, "rb");
 241     if (fp_remote == 0) {
 242         fprintf(stderr, "[download_from_remote] fail to open remote file %s\n",url);
 243         return NULL;
 244     }
 245     if ((fp = fopen(fn, "wb")) == 0) {
 246         fprintf(stderr, "[download_from_remote] fail to create file in the working directory %s\n",fn);
 247         knet_close(fp_remote);
 248         return NULL;
 249     }
 250     buf = (uint8_t*)calloc(buf_size, 1);
 251     while ((l = knet_read(fp_remote, buf, buf_size)) != 0)
 252         fwrite(buf, 1, l, fp);
 253     free(buf);
 254     fclose(fp);
 255     knet_close(fp_remote);
 256
 257     return fopen(fn, "r");
 258 }
 259 #endif
 260
 261 faidx_t *fai_load(const char *fn)
 262 {
 263         char *str;
 264         FILE *fp;
 265         faidx_t *fai;
 266         str = (char*)calloc(strlen(fn) + 5, 1);
 267         sprintf(str, "%s.fai", fn);
 268
 269 #ifdef _USE_KNETFILE
 270     if (strstr(fn, "ftp://") == fn || strstr(fn, "http://") == fn)
 271     {
 272         fp = download_and_open(str);
 273         if ( !fp )
 274         {
 275             fprintf(stderr, "[fai_load] failed to open remote FASTA index %s\n", str);
 276             free(str);
 277             return 0;
 278         }
 279     }
 280     else
 281 #endif
 282         fp = fopen(str, "rb");
 283         if (fp == 0) {
 284                 fprintf(stderr, "[fai_load] build FASTA index.\n");
 285                 fai_build(fn);
 286                 fp = fopen(str, "rb");
 287                 if (fp == 0) {
 288                         fprintf(stderr, "[fai_load] fail to open FASTA index.\n");
 289                         free(str);
 290                         return 0;
 291                 }
 292         }
 293
 294         fai = fai_read(fp);
 295         fclose(fp);
 296
 297         fai->rz = razf_open(fn, "rb");
 298         free(str);
 299         if (fai->rz == 0) {
 300                 fprintf(stderr, "[fai_load] fail to open FASTA file.\n");
 301                 return 0;
 302         }
 303         return fai;
 304 }
 305
 306 char *fai_fetch(const faidx_t *fai, const char *str, int *len)
 307 {
 308         char *s, *p, c;
 309         int i, l, k;
 310         khiter_t iter;
 311         faidx1_t val;
 312         khash_t(s) *h;
 313         int beg, end;
 314
 315         beg = end = -1;
 316         h = fai->hash;
 317         l = strlen(str);
 318         p = s = (char*)malloc(l+1);
 319         /* squeeze out "," */
 320         for (i = k = 0; i != l; ++i)
 321                 if (str[i] != ',' && !isspace(str[i])) s[k++] = str[i];
 322         s[k] = 0;
 323         for (i = 0; i != k; ++i) if (s[i] == ':') break;
 324         s[i] = 0;
 325         iter = kh_get(s, h, s); /* get the ref_id */
 326         if (iter == kh_end(h)) {
 327                 *len = 0;
 328                 free(s); return 0;
 329         }
 330         val = kh_value(h, iter);
 331         if (i == k) { /* dump the whole sequence */
 332                 beg = 0; end = val.len;
 333         } else {
 334                 for (p = s + i + 1; i != k; ++i) if (s[i] == '-') break;
 335                 beg = atoi(p);
 336                 if (i < k) {
 337                         p = s + i + 1;
 338                         end = atoi(p);
 339                 } else end = val.len;
 340         }
 341         if (beg > 0) --beg;
 342         if (beg >= val.len) beg = val.len;
 343         if (end >= val.len) end = val.len;
 344         if (beg > end) beg = end;
 345         free(s);
 346
 347         // now retrieve the sequence
 348         l = 0;
 349         s = (char*)malloc(end - beg + 2);
 350         razf_seek(fai->rz, val.offset + beg / val.line_blen * val.line_len + beg % val.line_blen, SEEK_SET);
 351         while (razf_read(fai->rz, &c, 1) == 1 && l < end - beg && !fai->rz->z_err)
 352                 if (isgraph(c)) s[l++] = c;
 353         s[l] = '\0';
 354         *len = l;
 355         return s;
 356 }
 357
 358 int faidx_main(int argc, char *argv[])
 359 {
 360         if (argc == 1) {
 361                 fprintf(stderr, "Usage: faidx <in.fasta> [<reg> [...]]\n");
 362                 return 1;
 363         } else {
 364                 if (argc == 2) fai_build(argv[1]);
 365                 else {
 366                         int i, j, k, l;
 367                         char *s;
 368                         faidx_t *fai;
 369                         fai = fai_load(argv[1]);
 370                         if (fai == 0) return 1;
 371                         for (i = 2; i != argc; ++i) {
 372                                 printf(">%s\n", argv[i]);
 373                                 s = fai_fetch(fai, argv[i], &l);
 374                                 for (j = 0; j < l; j += 60) {
 375                                         for (k = 0; k < 60 && k < l - j; ++k)
 376                                                 putchar(s[j + k]);
 377                                         putchar('\n');
 378                                 }
 379                                 free(s);
 380                         }
 381                         fai_destroy(fai);
 382                 }
 383         }
 384         return 0;
 385 }
 386
 387 int faidx_fetch_nseq(const faidx_t *fai)
 388 {
 389         return fai->n;
 390 }
 391
 392 char *faidx_fetch_seq(const faidx_t *fai, char *c_name, int p_beg_i, int p_end_i, int *len)
 393 {
 394         int l;
 395         char c;
 396     khiter_t iter;
 397     faidx1_t val;
 398         char *seq=NULL;
 399
 400     // Adjust position
 401     iter = kh_get(s, fai->hash, c_name);
 402     if(iter == kh_end(fai->hash)) return 0;
 403     val = kh_value(fai->hash, iter);
 404         if(p_end_i < p_beg_i) p_beg_i = p_end_i;
 405     if(p_beg_i < 0) p_beg_i = 0;
 406     else if(val.len <= p_beg_i) p_beg_i = val.len - 1;
 407     if(p_end_i < 0) p_end_i = 0;
 408     else if(val.len <= p_end_i) p_end_i = val.len - 1;
 409
 410     // Now retrieve the sequence
 411         l = 0;
 412         seq = (char*)malloc(p_end_i - p_beg_i + 2);
 413         razf_seek(fai->rz, val.offset + p_beg_i / val.line_blen * val.line_len + p_beg_i % val.line_blen, SEEK_SET);
 414         while (razf_read(fai->rz, &c, 1) == 1 && l < p_end_i - p_beg_i + 1)
 415                 if (isgraph(c)) seq[l++] = c;
 416         seq[l] = '\0';
 417         *len = l;
 418         return seq;
 419 }
 420
 421 #ifdef FAIDX_MAIN
 422 int main(int argc, char *argv[]) { return faidx_main(argc, argv); }
 423 #endif