samtools/faidx.c.pysam.c

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