annotate ezBAMQC/src/htslib/faidx.c @ 1:b4ff32b63fcd

Uploaded
author cshl-bsr
date Tue, 29 Mar 2016 15:25:51 -0400
parents dfa3745e5fd8
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
0
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
1 /* faidx.c -- FASTA random access.
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
2
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
3 Copyright (C) 2008, 2009, 2013-2015 Genome Research Ltd.
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
4 Portions copyright (C) 2011 Broad Institute.
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
5
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
6 Author: Heng Li <lh3@sanger.ac.uk>
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
7
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
8 Permission is hereby granted, free of charge, to any person obtaining a copy
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
9 of this software and associated documentation files (the "Software"), to deal
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
10 in the Software without restriction, including without limitation the rights
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
11 to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
12 copies of the Software, and to permit persons to whom the Software is
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
13 furnished to do so, subject to the following conditions:
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
14
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
15 The above copyright notice and this permission notice shall be included in
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
16 all copies or substantial portions of the Software.
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
17
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
18 THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
19 IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
20 FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
21 THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
22 LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
23 FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
24 DEALINGS IN THE SOFTWARE. */
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
25
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
26 #include <ctype.h>
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
27 #include <string.h>
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
28 #include <stdlib.h>
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
29 #include <stdio.h>
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
30 #include <stdint.h>
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
31
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
32 #include "htslib/bgzf.h"
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
33 #include "htslib/faidx.h"
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
34 #include "htslib/hfile.h"
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
35 #include "htslib/khash.h"
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
36
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
37 typedef struct {
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
38 int32_t line_len, line_blen;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
39 int64_t len;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
40 uint64_t offset;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
41 } faidx1_t;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
42 KHASH_MAP_INIT_STR(s, faidx1_t)
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
43
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
44 struct __faidx_t {
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
45 BGZF *bgzf;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
46 int n, m;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
47 char **name;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
48 khash_t(s) *hash;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
49 };
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
50
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
51 #ifndef kroundup32
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
52 #define kroundup32(x) (--(x), (x)|=(x)>>1, (x)|=(x)>>2, (x)|=(x)>>4, (x)|=(x)>>8, (x)|=(x)>>16, ++(x))
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
53 #endif
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
54
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
55 static inline void fai_insert_index(faidx_t *idx, const char *name, int len, int line_len, int line_blen, uint64_t offset)
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
56 {
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
57 khint_t k;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
58 int ret;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
59 faidx1_t t;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
60 if (idx->n == idx->m) {
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
61 idx->m = idx->m? idx->m<<1 : 16;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
62 idx->name = (char**)realloc(idx->name, sizeof(char*) * idx->m);
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
63 }
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
64 idx->name[idx->n] = strdup(name);
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
65 k = kh_put(s, idx->hash, idx->name[idx->n], &ret);
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
66 t.len = len; t.line_len = line_len; t.line_blen = line_blen; t.offset = offset;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
67 kh_value(idx->hash, k) = t;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
68 ++idx->n;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
69 }
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
70
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
71 faidx_t *fai_build_core(BGZF *bgzf)
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
72 {
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
73 char *name;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
74 int c;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
75 int l_name, m_name;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
76 int line_len, line_blen, state;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
77 int l1, l2;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
78 faidx_t *idx;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
79 uint64_t offset;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
80 int64_t len;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
81
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
82 idx = (faidx_t*)calloc(1, sizeof(faidx_t));
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
83 idx->hash = kh_init(s);
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
84 name = 0; l_name = m_name = 0;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
85 len = line_len = line_blen = -1; state = 0; l1 = l2 = -1; offset = 0;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
86 while ( (c=bgzf_getc(bgzf))>=0 ) {
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
87 if (c == '\n') { // an empty line
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
88 if (state == 1) {
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
89 offset = bgzf_utell(bgzf);
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
90 continue;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
91 } else if ((state == 0 && len < 0) || state == 2) continue;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
92 }
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
93 if (c == '>') { // fasta header
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
94 if (len >= 0)
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
95 fai_insert_index(idx, name, len, line_len, line_blen, offset);
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
96 l_name = 0;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
97 while ( (c=bgzf_getc(bgzf))>=0 && !isspace(c)) {
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
98 if (m_name < l_name + 2) {
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
99 m_name = l_name + 2;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
100 kroundup32(m_name);
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
101 name = (char*)realloc(name, m_name);
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
102 }
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
103 name[l_name++] = c;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
104 }
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
105 name[l_name] = '\0';
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
106 if ( c<0 ) {
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
107 fprintf(stderr, "[fai_build_core] the last entry has no sequence\n");
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
108 free(name); fai_destroy(idx);
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
109 return 0;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
110 }
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
111 if (c != '\n') while ( (c=bgzf_getc(bgzf))>=0 && c != '\n');
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
112 state = 1; len = 0;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
113 offset = bgzf_utell(bgzf);
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
114 } else {
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
115 if (state == 3) {
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
116 fprintf(stderr, "[fai_build_core] inlined empty line is not allowed in sequence '%s'.\n", name);
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
117 free(name); fai_destroy(idx);
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
118 return 0;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
119 }
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
120 if (state == 2) state = 3;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
121 l1 = l2 = 0;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
122 do {
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
123 ++l1;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
124 if (isgraph(c)) ++l2;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
125 } while ( (c=bgzf_getc(bgzf))>=0 && c != '\n');
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
126 if (state == 3 && l2) {
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
127 fprintf(stderr, "[fai_build_core] different line length in sequence '%s'.\n", name);
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
128 free(name); fai_destroy(idx);
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
129 return 0;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
130 }
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
131 ++l1; len += l2;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
132 if (state == 1) line_len = l1, line_blen = l2, state = 0;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
133 else if (state == 0) {
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
134 if (l1 != line_len || l2 != line_blen) state = 2;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
135 }
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
136 }
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
137 }
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
138 if ( name )
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
139 fai_insert_index(idx, name, len, line_len, line_blen, offset);
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
140 else
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
141 {
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
142 free(idx);
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
143 return NULL;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
144 }
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
145 free(name);
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
146 return idx;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
147 }
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
148
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
149 void fai_save(const faidx_t *fai, FILE *fp)
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
150 {
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
151 khint_t k;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
152 int i;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
153 for (i = 0; i < fai->n; ++i) {
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
154 faidx1_t x;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
155 k = kh_get(s, fai->hash, fai->name[i]);
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
156 x = kh_value(fai->hash, k);
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
157 #ifdef _WIN32
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
158 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);
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
159 #else
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
160 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);
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
161 #endif
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
162 }
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
163 }
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
164
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
165 faidx_t *fai_read(FILE *fp)
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
166 {
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
167 faidx_t *fai;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
168 char *buf, *p;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
169 int len, line_len, line_blen;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
170 #ifdef _WIN32
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
171 long offset;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
172 #else
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
173 long long offset;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
174 #endif
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
175 fai = (faidx_t*)calloc(1, sizeof(faidx_t));
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
176 fai->hash = kh_init(s);
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
177 buf = (char*)calloc(0x10000, 1);
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
178 while (!feof(fp) && fgets(buf, 0x10000, fp)) {
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
179 for (p = buf; *p && isgraph(*p); ++p);
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
180 *p = 0; ++p;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
181 #ifdef _WIN32
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
182 sscanf(p, "%d%ld%d%d", &len, &offset, &line_blen, &line_len);
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
183 #else
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
184 sscanf(p, "%d%lld%d%d", &len, &offset, &line_blen, &line_len);
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
185 #endif
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
186 fai_insert_index(fai, buf, len, line_len, line_blen, offset);
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
187 }
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
188 free(buf);
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
189 return fai;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
190 }
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
191
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
192 void fai_destroy(faidx_t *fai)
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
193 {
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
194 int i;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
195 for (i = 0; i < fai->n; ++i) free(fai->name[i]);
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
196 free(fai->name);
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
197 kh_destroy(s, fai->hash);
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
198 if (fai->bgzf) bgzf_close(fai->bgzf);
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
199 free(fai);
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
200 }
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
201
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
202 int fai_build(const char *fn)
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
203 {
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
204 char *str;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
205 BGZF *bgzf;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
206 FILE *fp;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
207 faidx_t *fai;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
208 str = (char*)calloc(strlen(fn) + 5, 1);
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
209 sprintf(str, "%s.fai", fn);
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
210 bgzf = bgzf_open(fn, "r");
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
211 if ( !bgzf ) {
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
212 fprintf(stderr, "[fai_build] fail to open the FASTA file %s\n",fn);
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
213 free(str);
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
214 return -1;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
215 }
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
216 if ( bgzf->is_compressed ) bgzf_index_build_init(bgzf);
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
217 fai = fai_build_core(bgzf);
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
218 if ( !fai )
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
219 {
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
220 if ( bgzf->is_compressed && bgzf->is_gzip ) fprintf(stderr,"Cannot index files compressed with gzip, please use bgzip\n");
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
221 free(str);
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
222 return -1;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
223 }
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
224 if ( bgzf->is_compressed ) bgzf_index_dump(bgzf, fn, ".gzi");
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
225 bgzf_close(bgzf);
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
226 fp = fopen(str, "wb");
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
227 if ( !fp ) {
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
228 fprintf(stderr, "[fai_build] fail to write FASTA index %s\n",str);
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
229 fai_destroy(fai); free(str);
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
230 return -1;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
231 }
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
232 fai_save(fai, fp);
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
233 fclose(fp);
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
234 free(str);
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
235 fai_destroy(fai);
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
236 return 0;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
237 }
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
238
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
239 static FILE *download_and_open(const char *fn)
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
240 {
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
241 const int buf_size = 1 * 1024 * 1024;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
242 uint8_t *buf;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
243 FILE *fp;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
244 hFILE *fp_remote;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
245 const char *url = fn;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
246 const char *p;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
247 int l = strlen(fn);
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
248 for (p = fn + l - 1; p >= fn; --p)
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
249 if (*p == '/') break;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
250 fn = p + 1;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
251
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
252 // First try to open a local copy
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
253 fp = fopen(fn, "r");
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
254 if (fp)
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
255 return fp;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
256
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
257 // If failed, download from remote and open
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
258 fp_remote = hopen(url, "rb");
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
259 if (fp_remote == 0) {
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
260 fprintf(stderr, "[download_from_remote] fail to open remote file %s\n",url);
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
261 return NULL;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
262 }
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
263 if ((fp = fopen(fn, "wb")) == 0) {
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
264 fprintf(stderr, "[download_from_remote] fail to create file in the working directory %s\n",fn);
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
265 hclose_abruptly(fp_remote);
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
266 return NULL;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
267 }
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
268 buf = (uint8_t*)calloc(buf_size, 1);
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
269 while ((l = hread(fp_remote, buf, buf_size)) > 0)
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
270 fwrite(buf, 1, l, fp);
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
271 free(buf);
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
272 fclose(fp);
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
273 if (hclose(fp_remote) != 0)
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
274 fprintf(stderr, "[download_from_remote] fail to close remote file %s\n", url);
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
275
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
276 return fopen(fn, "r");
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
277 }
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
278
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
279 faidx_t *fai_load(const char *fn)
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
280 {
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
281 char *str;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
282 FILE *fp;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
283 faidx_t *fai;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
284 str = (char*)calloc(strlen(fn) + 5, 1);
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
285 sprintf(str, "%s.fai", fn);
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
286
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
287 if (hisremote(str))
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
288 {
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
289 fp = download_and_open(str);
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
290 if ( !fp )
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
291 {
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
292 fprintf(stderr, "[fai_load] failed to open remote FASTA index %s\n", str);
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
293 free(str);
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
294 return 0;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
295 }
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
296 }
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
297 else
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
298 fp = fopen(str, "rb");
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
299
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
300 if (fp == 0) {
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
301 fprintf(stderr, "[fai_load] build FASTA index.\n");
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
302 fai_build(fn);
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
303 fp = fopen(str, "rb");
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
304 if (fp == 0) {
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
305 fprintf(stderr, "[fai_load] fail to open FASTA index.\n");
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
306 free(str);
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
307 return 0;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
308 }
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
309 }
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
310
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
311 fai = fai_read(fp);
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
312 fclose(fp);
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
313
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
314 fai->bgzf = bgzf_open(fn, "rb");
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
315 free(str);
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
316 if (fai->bgzf == 0) {
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
317 fprintf(stderr, "[fai_load] fail to open FASTA file.\n");
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
318 return 0;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
319 }
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
320 if ( fai->bgzf->is_compressed==1 )
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
321 {
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
322 if ( bgzf_index_load(fai->bgzf, fn, ".gzi") < 0 )
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
323 {
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
324 fprintf(stderr, "[fai_load] failed to load .gzi index: %s[.gzi]\n", fn);
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
325 fai_destroy(fai);
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
326 return NULL;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
327 }
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
328 }
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
329 return fai;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
330 }
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
331
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
332 char *fai_fetch(const faidx_t *fai, const char *str, int *len)
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
333 {
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
334 char *s;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
335 int c, i, l, k, name_end;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
336 khiter_t iter;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
337 faidx1_t val;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
338 khash_t(s) *h;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
339 int beg, end;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
340
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
341 beg = end = -1;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
342 h = fai->hash;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
343 name_end = l = strlen(str);
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
344 s = (char*)malloc(l+1);
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
345 // remove space
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
346 for (i = k = 0; i < l; ++i)
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
347 if (!isspace(str[i])) s[k++] = str[i];
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
348 s[k] = 0; l = k;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
349 // determine the sequence name
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
350 for (i = l - 1; i >= 0; --i) if (s[i] == ':') break; // look for colon from the end
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
351 if (i >= 0) name_end = i;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
352 if (name_end < l) { // check if this is really the end
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
353 int n_hyphen = 0;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
354 for (i = name_end + 1; i < l; ++i) {
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
355 if (s[i] == '-') ++n_hyphen;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
356 else if (!isdigit(s[i]) && s[i] != ',') break;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
357 }
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
358 if (i < l || n_hyphen > 1) name_end = l; // malformated region string; then take str as the name
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
359 s[name_end] = 0;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
360 iter = kh_get(s, h, s);
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
361 if (iter == kh_end(h)) { // cannot find the sequence name
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
362 iter = kh_get(s, h, str); // try str as the name
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
363 if (iter == kh_end(h)) {
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
364 *len = 0;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
365 free(s); return 0;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
366 } else s[name_end] = ':', name_end = l;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
367 }
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
368 } else iter = kh_get(s, h, str);
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
369 if(iter == kh_end(h)) {
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
370 fprintf(stderr, "[fai_fetch] Warning - Reference %s not found in FASTA file, returning empty sequence\n", str);
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
371 free(s);
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
372 *len = -2;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
373 return 0;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
374 };
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
375 val = kh_value(h, iter);
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
376 // parse the interval
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
377 if (name_end < l) {
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
378 for (i = k = name_end + 1; i < l; ++i)
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
379 if (s[i] != ',') s[k++] = s[i];
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
380 s[k] = 0;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
381 beg = atoi(s + name_end + 1);
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
382 for (i = name_end + 1; i != k; ++i) if (s[i] == '-') break;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
383 end = i < k? atoi(s + i + 1) : val.len;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
384 if (beg > 0) --beg;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
385 } else beg = 0, end = val.len;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
386 if (beg >= val.len) beg = val.len;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
387 if (end >= val.len) end = val.len;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
388 if (beg > end) beg = end;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
389 free(s);
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
390
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
391 // now retrieve the sequence
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
392 int ret = bgzf_useek(fai->bgzf, val.offset + beg / val.line_blen * val.line_len + beg % val.line_blen, SEEK_SET);
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
393 if ( ret<0 )
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
394 {
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
395 *len = -1;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
396 fprintf(stderr, "[fai_fetch] Error: fai_fetch failed. (Seeking in a compressed, .gzi unindexed, file?)\n");
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
397 return NULL;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
398 }
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
399 l = 0;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
400 s = (char*)malloc(end - beg + 2);
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
401 while ( (c=bgzf_getc(fai->bgzf))>=0 && l < end - beg )
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
402 if (isgraph(c)) s[l++] = c;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
403 s[l] = '\0';
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
404 *len = l;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
405 return s;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
406 }
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
407
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
408 int faidx_fetch_nseq(const faidx_t *fai)
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
409 {
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
410 return fai->n;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
411 }
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
412
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
413 int faidx_nseq(const faidx_t *fai)
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
414 {
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
415 return fai->n;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
416 }
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
417
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
418 const char *faidx_iseq(const faidx_t *fai, int i)
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
419 {
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
420 return fai->name[i];
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
421 }
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
422
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
423 int faidx_seq_len(const faidx_t *fai, const char *seq)
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
424 {
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
425 khint_t k = kh_get(s, fai->hash, seq);
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
426 if ( k == kh_end(fai->hash) ) return -1;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
427 return kh_val(fai->hash, k).len;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
428 }
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
429
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
430 char *faidx_fetch_seq(const faidx_t *fai, const char *c_name, int p_beg_i, int p_end_i, int *len)
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
431 {
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
432 int l, c;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
433 khiter_t iter;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
434 faidx1_t val;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
435 char *seq=NULL;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
436
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
437 // Adjust position
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
438 iter = kh_get(s, fai->hash, c_name);
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
439 if (iter == kh_end(fai->hash))
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
440 {
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
441 *len = -2;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
442 fprintf(stderr, "[fai_fetch_seq] The sequence \"%s\" not found\n", c_name);
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
443 return NULL;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
444 }
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
445 val = kh_value(fai->hash, iter);
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
446 if(p_end_i < p_beg_i) p_beg_i = p_end_i;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
447 if(p_beg_i < 0) p_beg_i = 0;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
448 else if(val.len <= p_beg_i) p_beg_i = val.len - 1;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
449 if(p_end_i < 0) p_end_i = 0;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
450 else if(val.len <= p_end_i) p_end_i = val.len - 1;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
451
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
452 // Now retrieve the sequence
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
453 int ret = bgzf_useek(fai->bgzf, val.offset + p_beg_i / val.line_blen * val.line_len + p_beg_i % val.line_blen, SEEK_SET);
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
454 if ( ret<0 )
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
455 {
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
456 *len = -1;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
457 fprintf(stderr, "[fai_fetch_seq] Error: fai_fetch failed. (Seeking in a compressed, .gzi unindexed, file?)\n");
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
458 return NULL;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
459 }
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
460 l = 0;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
461 seq = (char*)malloc(p_end_i - p_beg_i + 2);
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
462 while ( (c=bgzf_getc(fai->bgzf))>=0 && l < p_end_i - p_beg_i + 1)
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
463 if (isgraph(c)) seq[l++] = c;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
464 seq[l] = '\0';
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
465 *len = l;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
466 return seq;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
467 }
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
468
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
469 int faidx_has_seq(const faidx_t *fai, const char *seq)
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
470 {
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
471 khiter_t iter = kh_get(s, fai->hash, seq);
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
472 if (iter == kh_end(fai->hash)) return 0;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
473 return 1;
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
474 }
dfa3745e5fd8 Uploaded
youngkim
parents:
diff changeset
475