- char **list;
- int i, j, n, *fields, max_fields;
- khash_t(64) *hash;
- bam_init_header_hash(h);
- list = __bam_get_lines(fn, &n);
- hash = kh_init(64);
- max_fields = 0; fields = 0;
- for (i = 0; i < n; ++i) {
- char *str = list[i];
- int chr, n_fields, ret;
- khint_t k;
- uint64_t x;
- n_fields = ksplit_core(str, 0, &max_fields, &fields);
- if (n_fields < 2) continue;
- chr = bam_get_tid(h, str + fields[0]);
- if (chr < 0) {
- fprintf(stderr, "[load_pos] unknown reference sequence name: %s\n", str + fields[0]);
- continue;
- }
- x = (uint64_t)chr << 32 | (atoi(str + fields[1]) - 1);
- k = kh_put(64, hash, x, &ret);
- if (ret == 0) {
- fprintf(stderr, "[load_pos] position %s:%s has been loaded.\n", str+fields[0], str+fields[1]);
- continue;
- }
- kh_val(hash, k) = 0;
- if (n_fields > 2) {
- // count
- for (j = 2; j < n_fields; ++j) {
- char *s = str + fields[j];
- if ((*s != '+' && *s != '-') || !isdigit(s[1])) break;
- }
- if (j > 2) { // update kh_val()
- int *q, y, z;
- q = kh_val(hash, k) = (int*)calloc(j - 1, sizeof(int));
- q[0] = j - 2; z = j; y = 1;
- for (j = 2; j < z; ++j)
- q[y++] = atoi(str + fields[j]);
- }
- }
- free(str);
- }
- free(list); free(fields);
- return hash;
-}
-
-// an analogy to pileup_func() below
-static int glt3_func(uint32_t tid, uint32_t pos, int n, const bam_pileup1_t *pu, void *data)
-{
- pu_data_t *d = (pu_data_t*)data;
- bam_maqindel_ret_t *r = 0;
- int rb, *proposed_indels = 0;
- glf1_t *g;
- glf3_t *g3;
-
- if (d->fai == 0) {
- fprintf(stderr, "[glt3_func] reference sequence is required for generating GLT. Abort!\n");
- exit(1);
- }
- if (d->hash) { // only output a list of sites
- khint_t k = kh_get(64, d->hash, (uint64_t)tid<<32|pos);
- if (k == kh_end(d->hash)) return 0;
- proposed_indels = kh_val(d->hash, k);
- }
- g3 = glf3_init1();
- if (d->fai && (int)tid != d->tid) {
- if (d->ref) { // then write the end mark
- g3->rtype = GLF3_RTYPE_END;
- glf3_write1(d->fp_glf, g3);
- }
- glf3_ref_write(d->fp_glf, d->h->target_name[tid], d->h->target_len[tid]); // write reference
- free(d->ref);
- d->ref = fai_fetch(d->fai, d->h->target_name[tid], &d->len);
- d->tid = tid;
- d->last_pos = 0;
- }
- rb = (d->ref && (int)pos < d->len)? d->ref[pos] : 'N';
- g = bam_maqcns_glfgen(n, pu, bam_nt16_table[rb], d->c);
- memcpy(g3, g, sizeof(glf1_t));
- g3->rtype = GLF3_RTYPE_SUB;
- g3->offset = pos - d->last_pos;
- d->last_pos = pos;
- glf3_write1(d->fp_glf, g3);
- if (pos < d->len) {
- int m = (!d->max_depth || d->max_depth>n) ? n : d->max_depth;
- if (proposed_indels)
- r = bam_maqindel(m, pos, d->ido, pu, d->ref, proposed_indels[0], proposed_indels+1);
- else r = bam_maqindel(m, pos, d->ido, pu, d->ref, 0, 0);
- }
- if (r) { // then write indel line
- int het = 3 * n, min;
- min = het;
- if (min > r->gl[0]) min = r->gl[0];
- if (min > r->gl[1]) min = r->gl[1];
- g3->ref_base = 0;
- g3->rtype = GLF3_RTYPE_INDEL;
- memset(g3->lk, 0, 10);
- g3->lk[0] = r->gl[0] - min < 255? r->gl[0] - min : 255;
- g3->lk[1] = r->gl[1] - min < 255? r->gl[1] - min : 255;
- g3->lk[2] = het - min < 255? het - min : 255;
- g3->offset = 0;
- g3->indel_len[0] = r->indel1;
- g3->indel_len[1] = r->indel2;
- g3->min_lk = min < 255? min : 255;
- g3->max_len = (abs(r->indel1) > abs(r->indel2)? abs(r->indel1) : abs(r->indel2)) + 1;
- g3->indel_seq[0] = strdup(r->s[0]+1);
- g3->indel_seq[1] = strdup(r->s[1]+1);
- glf3_write1(d->fp_glf, g3);
- bam_maqindel_ret_destroy(r);
- }
- free(g);
- glf3_destroy1(g3);