From c350a570e955d5b4c2c13f607b7442df6d332c67 Mon Sep 17 00:00:00 2001
From: Heng Li <lh3@live.co.uk>
Date: Fri, 22 Apr 2011 13:46:42 +0000
Subject: [PATCH] Removed the pileup command

---
 Makefile        |  10 +-
 bam2bcf_indel.c |   5 +-
 bam_maqcns.c    | 626 ------------------------------------------------
 bam_maqcns.h    |  61 -----
 bam_plcmd.c     | 465 -----------------------------------
 bam_tview.c     |  34 ++-
 bamtk.c         |   6 -
 glf.c           | 236 ------------------
 glf.h           |  56 -----
 9 files changed, 32 insertions(+), 1467 deletions(-)
 delete mode 100644 bam_maqcns.c
 delete mode 100644 bam_maqcns.h
 delete mode 100644 glf.c
 delete mode 100644 glf.h

diff --git a/Makefile b/Makefile
index 4fc03d5..1b03921 100644
--- a/Makefile
+++ b/Makefile
@@ -3,9 +3,9 @@ CFLAGS=		-g -Wall -O2 #-m64 #-arch ppc
 DFLAGS=		-D_FILE_OFFSET_BITS=64 -D_LARGEFILE64_SOURCE -D_USE_KNETFILE -D_CURSES_LIB=1
 KNETFILE_O=	knetfile.o
 LOBJS=		bgzf.o kstring.o bam_aux.o bam.o bam_import.o sam.o bam_index.o	\
-			bam_pileup.o bam_lpileup.o bam_md.o glf.o razf.o faidx.o bedidx.o \
+			bam_pileup.o bam_lpileup.o bam_md.o razf.o faidx.o bedidx.o \
 			$(KNETFILE_O) bam_sort.o sam_header.o bam_reheader.o kprobaln.o bam_cat.o
-AOBJS=		bam_tview.o bam_maqcns.o bam_plcmd.o sam_view.o	\
+AOBJS=		bam_tview.o bam_plcmd.o sam_view.o	\
 			bam_rmdup.o bam_rmdupse.o bam_mate.o bam_stat.o bam_color.o	\
 			bamtk.o kaln.o bam2bcf.o bam2bcf_indel.o errmod.o sample.o \
 			cut_target.o phase.o bam2depth.o
@@ -54,14 +54,12 @@ bam.o:bam.h razf.h bam_endian.h kstring.h sam_header.h
 sam.o:sam.h bam.h
 bam_import.o:bam.h kseq.h khash.h razf.h
 bam_pileup.o:bam.h razf.h ksort.h
-bam_plcmd.o:bam.h faidx.h bam_maqcns.h glf.h bcftools/bcf.h bam2bcf.h
+bam_plcmd.o:bam.h faidx.h bcftools/bcf.h bam2bcf.h
 bam_index.o:bam.h khash.h ksort.h razf.h bam_endian.h
 bam_lpileup.o:bam.h ksort.h
-bam_tview.o:bam.h faidx.h bam_maqcns.h
-bam_maqcns.o:bam.h ksort.h bam_maqcns.h kaln.h
+bam_tview.o:bam.h faidx.h
 bam_sort.o:bam.h ksort.h razf.h
 bam_md.o:bam.h faidx.h
-glf.o:glf.h
 sam_header.o:sam_header.h khash.h
 bcf.o:bcftools/bcf.h
 bam2bcf.o:bam2bcf.h errmod.h bcftools/bcf.h
diff --git a/bam2bcf_indel.c b/bam2bcf_indel.c
index 899428f..f37a783 100644
--- a/bam2bcf_indel.c
+++ b/bam2bcf_indel.c
@@ -3,12 +3,14 @@
 #include <string.h>
 #include "bam.h"
 #include "bam2bcf.h"
-#include "ksort.h"
 #include "kaln.h"
 #include "kprobaln.h"
 #include "khash.h"
 KHASH_SET_INIT_STR(rg)
 
+#include "ksort.h"
+KSORT_INIT_GENERIC(uint32_t)
+
 #define MINUS_CONST 0x10000000
 #define INDEL_WINDOW_SIZE 50
 
@@ -110,7 +112,6 @@ static inline int est_indelreg(int pos, const char *ref, int l, char *ins4)
 int bcf_call_gap_prep(int n, int *n_plp, bam_pileup1_t **plp, int pos, bcf_callaux_t *bca, const char *ref,
 					  const void *rghash)
 {
-	extern void ks_introsort_uint32_t(int, uint32_t*);
 	int i, s, j, k, t, n_types, *types, max_rd_len, left, right, max_ins, *score1, *score2, max_ref2;
 	int N, K, l_run, ref_type, n_alt;
 	char *inscns = 0, *ref2, *query, **ref_sample;
diff --git a/bam_maqcns.c b/bam_maqcns.c
deleted file mode 100644
index 9931036..0000000
--- a/bam_maqcns.c
+++ /dev/null
@@ -1,626 +0,0 @@
-#include <math.h>
-#include <assert.h>
-#include "bam.h"
-#include "bam_maqcns.h"
-#include "ksort.h"
-#include "errmod.h"
-#include "kaln.h"
-KSORT_INIT_GENERIC(uint32_t)
-
-#define INDEL_WINDOW_SIZE 50
-#define INDEL_EXT_DEP 0.9
-
-typedef struct __bmc_aux_t {
-	int max;
-	uint32_t *info;
-	uint16_t *info16;
-	errmod_t *em;
-} bmc_aux_t;
-
-typedef struct {
-	float esum[4], fsum[4];
-	uint32_t c[4];
-} glf_call_aux_t;
-
-/*
-  P(<b1,b2>) = \theta \sum_{i=1}^{N-1} 1/i
-  P(D|<b1,b2>) = \sum_{k=1}^{N-1} p_k 1/2 [(k/N)^n_2(1-k/N)^n_1 + (k/N)^n1(1-k/N)^n_2]
-  p_k = 1/k / \sum_{i=1}^{N-1} 1/i
- */
-static void cal_het(bam_maqcns_t *aa)
-{
-	int k, n1, n2;
-	double sum_harmo; // harmonic sum
-	double poly_rate;
-
-	free(aa->lhet);
-	aa->lhet = (double*)calloc(256 * 256, sizeof(double));
-	sum_harmo = 0.0;
-	for (k = 1; k <= aa->n_hap - 1; ++k)
-		sum_harmo += 1.0 / k;
-	for (n1 = 0; n1 < 256; ++n1) {
-		for (n2 = 0; n2 < 256; ++n2) {
-			long double sum = 0.0;
-			double lC = aa->errmod == BAM_ERRMOD_SOAP? 0 : lgamma(n1+n2+1) - lgamma(n1+1) - lgamma(n2+1);
-			for (k = 1; k <= aa->n_hap - 1; ++k) {
-				double pk = 1.0 / k / sum_harmo;
-				double log1 = log((double)k/aa->n_hap);
-				double log2 = log(1.0 - (double)k/aa->n_hap);
-				sum += pk * 0.5 * (expl(log1*n2) * expl(log2*n1) + expl(log1*n1) * expl(log2*n2));
-			}
-			aa->lhet[n1<<8|n2] = lC + logl(sum);
-		}
-	}
-	poly_rate = aa->het_rate * sum_harmo;
-	aa->q_r = -4.343 * log(2.0 * poly_rate / (1.0 - poly_rate));
-}
-
-/** initialize the helper structure */
-static void cal_coef(bam_maqcns_t *aa)
-{
-	int k, n, q;
-	long double sum_a[257], b[256], q_c[256], tmp[256], fk2[256];
-	double *lC;
-
-	if (aa->errmod == BAM_ERRMOD_MAQ2) return; // no need to do the following
-	// aa->lhet will be allocated and initialized 
-	free(aa->fk); free(aa->coef);
-	aa->coef = 0;
-	aa->fk = (double*)calloc(256, sizeof(double));
-	aa->fk[0] = fk2[0] = 1.0;
-	for (n = 1; n != 256; ++n) {
-		aa->fk[n] = pow(aa->theta, n) * (1.0 - aa->eta) + aa->eta;
-		fk2[n] = aa->fk[n>>1]; // this is an approximation, assuming reads equally likely come from both strands
-	}
-	if (aa->errmod == BAM_ERRMOD_SOAP) return;
-	aa->coef = (double*)calloc(256*256*64, sizeof(double));
-	lC = (double*)calloc(256 * 256, sizeof(double));
-	for (n = 1; n != 256; ++n)
-		for (k = 1; k <= n; ++k)
-			lC[n<<8|k] = lgamma(n+1) - lgamma(k+1) - lgamma(n-k+1);
-	for (q = 1; q != 64; ++q) {
-		double e = pow(10.0, -q/10.0);
-		double le = log(e);
-		double le1 = log(1.0-e);
-		for (n = 1; n != 256; ++n) {
-			double *coef = aa->coef + (q<<16|n<<8);
-			sum_a[n+1] = 0.0;
-			for (k = n; k >= 0; --k) { // a_k = \sum_{i=k}^n C^n_k \epsilon^k (1-\epsilon)^{n-k}
-				sum_a[k] = sum_a[k+1] + expl(lC[n<<8|k] + k*le + (n-k)*le1);
-				b[k] = sum_a[k+1] / sum_a[k];
-				if (b[k] > 0.99) b[k] = 0.99;
-			}
-			for (k = 0; k != n; ++k) // log(\bar\beta_{nk}(\bar\epsilon)^{f_k})
-				q_c[k] = -4.343 * fk2[k] * logl(b[k] / e);
-			for (k = 1; k != n; ++k) q_c[k] += q_c[k-1]; // \prod_{i=0}^k c_i
-			for (k = 0; k <= n; ++k) { // powl() in 64-bit mode seems broken on my Mac OS X 10.4.9
-				tmp[k] = -4.343 * logl(1.0 - expl(fk2[k] * logl(b[k])));
-				coef[k] = (k? q_c[k-1] : 0) + tmp[k]; // this is the final c_{nk}
-			}
-		}
-	}
-	free(lC);
-}
-
-bam_maqcns_t *bam_maqcns_init()
-{
-	bam_maqcns_t *bm;
-	bm = (bam_maqcns_t*)calloc(1, sizeof(bam_maqcns_t));
-	bm->aux = (bmc_aux_t*)calloc(1, sizeof(bmc_aux_t));
-	bm->het_rate = 0.001;
-	bm->theta = 0.83f;
-	bm->n_hap = 2;
-	bm->eta = 0.03;
-	bm->cap_mapQ = 60;
-	bm->min_baseQ = 13;
-	return bm;
-}
-
-void bam_maqcns_prepare(bam_maqcns_t *bm)
-{
-	if (bm->errmod == BAM_ERRMOD_MAQ2) bm->aux->em = errmod_init(1. - bm->theta);
-	cal_coef(bm); cal_het(bm);
-}
-
-void bam_maqcns_destroy(bam_maqcns_t *bm)
-{
-	if (bm == 0) return;
-	free(bm->lhet); free(bm->fk); free(bm->coef); free(bm->aux->info); free(bm->aux->info16);
-	if (bm->aux->em) errmod_destroy(bm->aux->em);
-	free(bm->aux); free(bm);
-}
-
-glf1_t *bam_maqcns_glfgen(int _n, const bam_pileup1_t *pl, uint8_t ref_base, bam_maqcns_t *bm)
-{
-	glf_call_aux_t *b = 0;
-	int i, j, k, w[8], c, n;
-	glf1_t *g = (glf1_t*)calloc(1, sizeof(glf1_t));
-	float p[16], min_p = 1e30;
-	uint64_t rms;
-
-	g->ref_base = ref_base;
-	if (_n == 0) return g;
-
-	// construct aux array
-	if (bm->aux->max < _n) {
-		bm->aux->max = _n;
-		kroundup32(bm->aux->max);
-		bm->aux->info = (uint32_t*)realloc(bm->aux->info, 4 * bm->aux->max);
-		bm->aux->info16 = (uint16_t*)realloc(bm->aux->info16, 2 * bm->aux->max);
-	}
-	for (i = n = 0, rms = 0; i < _n; ++i) {
-		const bam_pileup1_t *p = pl + i;
-		uint32_t q, x = 0, qq;
-		uint16_t y = 0;
-		if (p->is_del || p->is_refskip || (p->b->core.flag&BAM_FUNMAP)) continue;
-		q = (uint32_t)bam1_qual(p->b)[p->qpos];
-		if (q < bm->min_baseQ) continue;
-		x |= (uint32_t)bam1_strand(p->b) << 18 | q << 8 | p->b->core.qual;
-		y |= bam1_strand(p->b)<<4;
-		if (p->b->core.qual < q) q = p->b->core.qual;
-		c = p->b->core.qual < bm->cap_mapQ? p->b->core.qual : bm->cap_mapQ;
-		rms += c * c;
-		x |= q << 24;
-		y |= q << 5;
-		qq = bam1_seqi(bam1_seq(p->b), p->qpos);
-		q = bam_nt16_nt4_table[qq? qq : ref_base];
-		if (!p->is_del && !p->is_refskip && q < 4) x |= 1 << 21 | q << 16, y |= q;
-		bm->aux->info16[n] = y;
-		bm->aux->info[n++] = x;
-	}
-	rms = (uint8_t)(sqrt((double)rms / n) + .499);
-	if (bm->errmod == BAM_ERRMOD_MAQ2) {
-		errmod_cal(bm->aux->em, n, 4, bm->aux->info16, p);
-		goto goto_glf;
-	}
-	ks_introsort(uint32_t, n, bm->aux->info);
-	// generate esum and fsum
-	b = (glf_call_aux_t*)calloc(1, sizeof(glf_call_aux_t));
-	for (k = 0; k != 8; ++k) w[k] = 0;
-	for (j = n - 1; j >= 0; --j) { // calculate esum and fsum
-		uint32_t info = bm->aux->info[j];
-		if (info>>24 < 4 && (info>>8&0x3f) != 0) info = 4<<24 | (info&0xffffff);
-		k = info>>16&7;
-		if (info>>24 > 0) {
-			b->esum[k&3] += bm->fk[w[k]] * (info>>24);
-			b->fsum[k&3] += bm->fk[w[k]];
-			if (w[k] < 0xff) ++w[k];
-			++b->c[k&3];
-		}
-	}
-	// rescale ->c[]
-	for (j = c = 0; j != 4; ++j) c += b->c[j];
-	if (c > 255) {
-		for (j = 0; j != 4; ++j) b->c[j] = (int)(254.0 * b->c[j] / c + 0.5);
-		for (j = c = 0; j != 4; ++j) c += b->c[j];
-	}
-	if (bm->errmod == BAM_ERRMOD_MAQ) {
-		// generate likelihood
-		for (j = 0; j != 4; ++j) {
-			// homozygous
-			float tmp1, tmp3;
-			int tmp2, bar_e;
-			for (k = 0, tmp1 = tmp3 = 0.0, tmp2 = 0; k != 4; ++k) {
-				if (j == k) continue;
-				tmp1 += b->esum[k]; tmp2 += b->c[k]; tmp3 += b->fsum[k];
-			}
-			if (tmp2) {
-				bar_e = (int)(tmp1 / tmp3 + 0.5);
-				if (bar_e < 4) bar_e = 4; // should not happen
-				if (bar_e > 63) bar_e = 63;
-				p[j<<2|j] = tmp1 + bm->coef[bar_e<<16|c<<8|tmp2];
-			} else p[j<<2|j] = 0.0; // all the bases are j
-			// heterozygous
-			for (k = j + 1; k < 4; ++k) {
-				for (i = 0, tmp2 = 0, tmp1 = tmp3 = 0.0; i != 4; ++i) {
-					if (i == j || i == k) continue;
-					tmp1 += b->esum[i]; tmp2 += b->c[i]; tmp3 += b->fsum[i];
-				}
-				if (tmp2) {
-					bar_e = (int)(tmp1 / tmp3 + 0.5);
-					if (bar_e < 4) bar_e = 4;
-					if (bar_e > 63) bar_e = 63;
-					p[j<<2|k] = p[k<<2|j] = -4.343 * bm->lhet[b->c[j]<<8|b->c[k]] + tmp1 + bm->coef[bar_e<<16|c<<8|tmp2];
-				} else p[j<<2|k] = p[k<<2|j] = -4.343 * bm->lhet[b->c[j]<<8|b->c[k]]; // all the bases are either j or k
-			}
-			//
-			for (k = 0; k != 4; ++k)
-				if (p[j<<2|k] < 0.0) p[j<<2|k] = 0.0;
-		}
-
-		{ // fix p[k<<2|k]
-			float max1, max2, min1, min2;
-			int max_k, min_k;
-			max_k = min_k = -1;
-			max1 = max2 = -1.0; min1 = min2 = 1e30;
-			for (k = 0; k < 4; ++k) {
-				if (b->esum[k] > max1) {
-					max2 = max1; max1 = b->esum[k]; max_k = k;
-				} else if (b->esum[k] > max2) max2 = b->esum[k];
-			}
-			for (k = 0; k < 4; ++k) {
-				if (p[k<<2|k] < min1) {
-					min2 = min1; min1 = p[k<<2|k]; min_k = k;
-				} else if (p[k<<2|k] < min2) min2 = p[k<<2|k];
-			}
-			if (max1 > max2 && (min_k != max_k || min1 + 1.0 > min2))
-				p[max_k<<2|max_k] = min1 > 1.0? min1 - 1.0 : 0.0;
-		}
-	} else if (bm->errmod == BAM_ERRMOD_SOAP) { // apply the SOAP model
-		// generate likelihood
-		for (j = 0; j != 4; ++j) {
-			float tmp;
-			// homozygous
-			for (k = 0, tmp = 0.0; k != 4; ++k)
-				if (j != k) tmp += b->esum[k];
-			p[j<<2|j] = tmp;
-			// heterozygous
-			for (k = j + 1; k < 4; ++k) {
-				for (i = 0, tmp = 0.0; i != 4; ++i)
-					if (i != j && i != k) tmp += b->esum[i];
-				p[j<<2|k] = p[k<<2|j] = -4.343 * bm->lhet[b->c[j]<<8|b->c[k]] + tmp;
-			}
-		}
-	}
-
-goto_glf:
-	// convert necessary information to glf1_t
-	g->ref_base = ref_base; g->max_mapQ = rms;
-	g->depth = n > 16777215? 16777215 : n;
-	for (j = 0; j != 4; ++j)
-		for (k = j; k < 4; ++k)
-			if (p[j<<2|k] < min_p) min_p = p[j<<2|k];
-	g->min_lk = min_p > 255.0? 255 : (int)(min_p + 0.5);
-	for (j = c = 0; j != 4; ++j)
-		for (k = j; k < 4; ++k)
-			g->lk[c++] = p[j<<2|k]-min_p > 255.0? 255 : (int)(p[j<<2|k]-min_p + 0.5);
-
-	free(b);
-	return g;
-}
-
-uint32_t glf2cns(const glf1_t *g, int q_r)
-{
-	int i, j, k, p[10], ref4;
-	uint32_t x = 0;
-	ref4 = bam_nt16_nt4_table[g->ref_base];
-	for (i = k = 0; i < 4; ++i)
-		for (j = i; j < 4; ++j) {
-			int prior = (i == ref4 && j == ref4? 0 : i == ref4 || j == ref4? q_r : q_r + 3);
-			p[k] = (g->lk[k] + prior)<<4 | i<<2 | j;
-			++k;
-		}
-	for (i = 1; i < 10; ++i) // insertion sort
-		for (j = i; j > 0 && p[j] < p[j-1]; --j)
-			k = p[j], p[j] = p[j-1], p[j-1] = k;
-	x = (1u<<(p[0]&3) | 1u<<(p[0]>>2&3)) << 28; // the best genotype
-	x |= (uint32_t)g->max_mapQ << 16; // rms mapQ
-	x |= ((p[1]>>4) - (p[0]>>4) < 256? (p[1]>>4) - (p[0]>>4) : 255) << 8; // consensus Q
-	for (k = 0; k < 10; ++k)
-		if ((p[k]&0xf) == (ref4<<2|ref4)) break;
-	if (k == 10) k = 9;
-	x |= (p[k]>>4) - (p[0]>>4) < 256? (p[k]>>4) - (p[0]>>4) : 255; // snp Q
-	return x;
-}
-
-uint32_t bam_maqcns_call(int n, const bam_pileup1_t *pl, bam_maqcns_t *bm)
-{
-	glf1_t *g;
-	uint32_t x;
-	if (n) {
-		g = bam_maqcns_glfgen(n, pl, 0xf, bm);
-		x = g->depth == 0? (0xfU<<28 | 0xfU<<24) : glf2cns(g, (int)(bm->q_r + 0.5));
-		free(g);
-	} else x = 0xfU<<28 | 0xfU<<24;
-	return x;
-}
-
-/************** *****************/
-
-bam_maqindel_opt_t *bam_maqindel_opt_init()
-{
-	bam_maqindel_opt_t *mi = (bam_maqindel_opt_t*)calloc(1, sizeof(bam_maqindel_opt_t));
-	mi->q_indel = 40;
-	mi->r_indel = 0.00015;
-	mi->r_snp = 0.001;
-	//
-	mi->mm_penalty = 3;
-	mi->indel_err = 4;
-	mi->ambi_thres = 10;
-	return mi;
-}
-
-void bam_maqindel_ret_destroy(bam_maqindel_ret_t *mir)
-{
-	if (mir == 0) return;
-	free(mir->s[0]); free(mir->s[1]); free(mir);
-}
-
-int bam_tpos2qpos(const bam1_core_t *c, const uint32_t *cigar, int32_t tpos, int is_left, int32_t *_tpos)
-{
-	int k, x = c->pos, y = 0, last_y = 0;
-	*_tpos = c->pos;
-	for (k = 0; k < c->n_cigar; ++k) {
-		int op = cigar[k] & BAM_CIGAR_MASK;
-		int l = cigar[k] >> BAM_CIGAR_SHIFT;
-		if (op == BAM_CMATCH) {
-			if (c->pos > tpos) return y;
-			if (x + l > tpos) {
-				*_tpos = tpos;
-				return y + (tpos - x);
-			}
-			x += l; y += l;
-			last_y = y;
-		} else if (op == BAM_CINS || op == BAM_CSOFT_CLIP) y += l;
-		else if (op == BAM_CDEL || op == BAM_CREF_SKIP) {
-			if (x + l > tpos) {
-				*_tpos = is_left? x : x + l;
-				return y;
-			}
-			x += l;
-		}
-	}
-	*_tpos = x;
-	return last_y;
-}
-
-#define MINUS_CONST 0x10000000
-
-bam_maqindel_ret_t *bam_maqindel(int n, int pos, const bam_maqindel_opt_t *mi, const bam_pileup1_t *pl, const char *ref,
-								 int _n_types, int *_types)
-{
-	int i, j, n_types, *types, left, right, max_rd_len = 0;
-	bam_maqindel_ret_t *ret = 0;
-	// if there is no proposed indel, check if there is an indel from the alignment
-	if (_n_types == 0) {
-		for (i = 0; i < n; ++i) {
-			const bam_pileup1_t *p = pl + i;
-			if (!(p->b->core.flag&BAM_FUNMAP) && p->indel != 0) break;
-		}
-		if (i == n) return 0; // no indel
-	}
-	{ // calculate how many types of indels are available (set n_types and types)
-		int m;
-		uint32_t *aux;
-		aux = (uint32_t*)calloc(n + _n_types + 1, 4);
-		m = 0;
-		aux[m++] = MINUS_CONST; // zero indel is always a type
-		for (i = 0; i < n; ++i) {
-			const bam_pileup1_t *p = pl + i;
-			if (!(p->b->core.flag&BAM_FUNMAP) && p->indel != 0)
-				aux[m++] = MINUS_CONST + p->indel;
-			j = bam_cigar2qlen(&p->b->core, bam1_cigar(p->b));
-			if (j > max_rd_len) max_rd_len = j;
-		}
-		if (_n_types) // then also add this to aux[]
-			for (i = 0; i < _n_types; ++i)
-				if (_types[i]) aux[m++] = MINUS_CONST + _types[i];
-		ks_introsort(uint32_t, m, aux);
-		// squeeze out identical types
-		for (i = 1, n_types = 1; i < m; ++i)
-			if (aux[i] != aux[i-1]) ++n_types;
-		types = (int*)calloc(n_types, sizeof(int));
-		j = 0;
-		types[j++] = aux[0] - MINUS_CONST; 
-		for (i = 1; i < m; ++i) {
-			if (aux[i] != aux[i-1])
-				types[j++] = aux[i] - MINUS_CONST;
-		}
-		free(aux);
-	}
-	{ // calculate left and right boundary
-		left = pos > INDEL_WINDOW_SIZE? pos - INDEL_WINDOW_SIZE : 0;
-		right = pos + INDEL_WINDOW_SIZE;
-		if (types[0] < 0) right -= types[0];
-		// in case the alignments stand out the reference
-		for (i = pos; i < right; ++i)
-			if (ref[i] == 0) break;
-		right = i;
-	}
-	{ // the core part
-		char *ref2, *rs, *inscns = 0;
-		int qr_snp, k, l, *score, *pscore, max_ins = types[n_types-1];
-		qr_snp = (int)(-4.343 * log(mi->r_snp) + .499);
-		if (max_ins > 0) { // get the consensus of inserted sequences
-			int *inscns_aux = (int*)calloc(4 * n_types * max_ins, sizeof(int));
-			// count occurrences
-			for (i = 0; i < n_types; ++i) {
-				if (types[i] <= 0) continue; // not insertion
-				for (j = 0; j < n; ++j) {
-					const bam_pileup1_t *p = pl + j;
-					if (!(p->b->core.flag&BAM_FUNMAP) && p->indel == types[i]) {
-						for (k = 1; k <= p->indel; ++k) {
-							int c = bam_nt16_nt4_table[bam1_seqi(bam1_seq(p->b), p->qpos + k)];
-							if (c < 4) ++inscns_aux[i*max_ins*4 + (k-1)*4 + c];
-						}
-					}
-				}
-			}
-			// construct the consensus of inserted sequence
-			inscns = (char*)calloc(n_types * max_ins, sizeof(char));
-			for (i = 0; i < n_types; ++i) {
-				for (j = 0; j < types[i]; ++j) {
-					int max = 0, max_k = -1, *ia = inscns_aux + i*max_ins*4 + j*4;
-					for (k = 0; k < 4; ++k) {
-						if (ia[k] > max) {
-							max = ia[k];
-							max_k = k;
-						}
-					}
-					inscns[i*max_ins + j] = max? 1<<max_k : 15;
-				}
-			}
-			free(inscns_aux);
-		}
-		// calculate score
-		ref2 = (char*)calloc(right - left + types[n_types-1] + 2, 1);
-		rs   = (char*)calloc(right - left + max_rd_len + types[n_types-1] + 2, 1);
-		score = (int*)calloc(n_types * n, sizeof(int));
-		pscore = (int*)calloc(n_types * n, sizeof(int));
-		for (i = 0; i < n_types; ++i) {
-			ka_param_t ap = ka_param_blast;
-			ap.band_width = 2 * types[n_types - 1] + 2;
-			ap.gap_end_ext = 0;
-			// write ref2
-			for (k = 0, j = left; j <= pos; ++j)
-				ref2[k++] = bam_nt16_nt4_table[bam_nt16_table[(int)ref[j]]];
-			if (types[i] <= 0) j += -types[i];
-			else for (l = 0; l < types[i]; ++l)
-					 ref2[k++] = bam_nt16_nt4_table[(int)inscns[i*max_ins + l]];
-			if (types[0] < 0) { // mask deleted sequences
-				int jj, tmp = types[i] >= 0? -types[0] : -types[0] + types[i];
-				for (jj = 0; jj < tmp && j < right && ref[j]; ++jj, ++j)
-					ref2[k++] = 4;
-			}
-			for (; j < right && ref[j]; ++j)
-				ref2[k++] = bam_nt16_nt4_table[bam_nt16_table[(int)ref[j]]];
-			if (j < right) right = j;
-			// calculate score for each read
-			for (j = 0; j < n; ++j) {
-				const bam_pileup1_t *p = pl + j;
-				int qbeg, qend, tbeg, tend;
-				if (p->b->core.flag & BAM_FUNMAP) continue;
-				qbeg = bam_tpos2qpos(&p->b->core, bam1_cigar(p->b), left,  0, &tbeg);
-				qend = bam_tpos2qpos(&p->b->core, bam1_cigar(p->b), right, 1, &tend);
-				assert(tbeg >= left);
-				for (l = qbeg; l < qend; ++l)
-					rs[l - qbeg] = bam_nt16_nt4_table[bam1_seqi(bam1_seq(p->b), l)];
-				{
-					int x, y, n_acigar, ps;
-					uint32_t *acigar;
-					ps = 0;
-					if (tend - tbeg + types[i] <= 0) {
-						score[i*n+j] = -(1<<20);
-						pscore[i*n+j] = 1<<20;
-						continue;
-					}
-					acigar = ka_global_core((uint8_t*)ref2 + tbeg - left, tend - tbeg + types[i], (uint8_t*)rs, qend - qbeg, &ap, &score[i*n+j], &n_acigar);
-					x = tbeg - left; y = 0;
-					for (l = 0; l < n_acigar; ++l) {
-						int op = acigar[l]&0xf;
-						int len = acigar[l]>>4;
-						if (op == BAM_CMATCH) {
-							int k;
-							for (k = 0; k < len; ++k)
-								if (ref2[x+k] != rs[y+k] && ref2[x+k] < 4)
-									ps += bam1_qual(p->b)[y+k] < qr_snp? bam1_qual(p->b)[y+k] : qr_snp;
-							x += len; y += len;
-						} else if (op == BAM_CINS || op == BAM_CSOFT_CLIP) {
-							if (op == BAM_CINS && l > 0 && l < n_acigar - 1) ps += mi->q_indel * len;
-							y += len;
-						} else if (op == BAM_CDEL) {
-							if (l > 0 && l < n_acigar - 1) ps += mi->q_indel * len;
-							x += len;
-						}
-					}
-					pscore[i*n+j] = ps;
-					/*if (1) { // for debugging only
-						fprintf(stderr, "id=%d, pos=%d, type=%d, j=%d, score=%d, psore=%d, %d, %d, %d, %d, %d, ",
-								j, pos+1, types[i], j, score[i*n+j], pscore[i*n+j], tbeg, tend, qbeg, qend, mi->q_indel);
-						for (l = 0; l < n_acigar; ++l) fprintf(stderr, "%d%c", acigar[l]>>4, "MIDS"[acigar[l]&0xf]);
-						fprintf(stderr, "\n");
-						for (l = 0; l < tend - tbeg + types[i]; ++l) fputc("ACGTN"[ref2[l+tbeg-left]], stderr);
-						fputc('\n', stderr);
-						for (l = 0; l < qend - qbeg; ++l) fputc("ACGTN"[rs[l]], stderr);
-						fputc('\n', stderr);
-						}*/
-					free(acigar);
-				}
-			}
-		}
-		{ // get final result
-			int *sum, max1, max2, max1_i, max2_i;
-			// pick up the best two score
-			sum = (int*)calloc(n_types, sizeof(int));
-			for (i = 0; i < n_types; ++i)
-				for (j = 0; j < n; ++j)
-					sum[i] += -pscore[i*n+j];
-			max1 = max2 = -0x7fffffff; max1_i = max2_i = -1;
-			for (i = 0; i < n_types; ++i) {
-				if (sum[i] > max1) {
-					max2 = max1; max2_i = max1_i; max1 = sum[i]; max1_i = i;
-				} else if (sum[i] > max2) {
-					max2 = sum[i]; max2_i = i;
-				}
-			}
-			free(sum);
-			// write ret
-			ret = (bam_maqindel_ret_t*)calloc(1, sizeof(bam_maqindel_ret_t));
-			ret->indel1 = types[max1_i]; ret->indel2 = types[max2_i];
-			ret->s[0] = (char*)calloc(abs(ret->indel1) + 2, 1);
-			ret->s[1] = (char*)calloc(abs(ret->indel2) + 2, 1);
-			// write indel sequence
-			if (ret->indel1 > 0) {
-				ret->s[0][0] = '+';
-				for (k = 0; k < ret->indel1; ++k)
-					ret->s[0][k+1] = bam_nt16_rev_table[(int)inscns[max1_i*max_ins + k]];
-			} else if (ret->indel1 < 0) {
-				ret->s[0][0] = '-';
-				for (k = 0; k < -ret->indel1 && ref[pos + k + 1]; ++k)
-					ret->s[0][k+1] = ref[pos + k + 1];
-			} else ret->s[0][0] = '*';
-			if (ret->indel2 > 0) {
-				ret->s[1][0] = '+';
-				for (k = 0; k < ret->indel2; ++k)
-					ret->s[1][k+1] = bam_nt16_rev_table[(int)inscns[max2_i*max_ins + k]];
-			} else if (ret->indel2 < 0) {
-				ret->s[1][0] = '-';
-				for (k = 0; k < -ret->indel2 && ref[pos + k + 1]; ++k)
-					ret->s[1][k+1] = ref[pos + k + 1];
-			} else ret->s[1][0] = '*';
-			// write count
-			for (i = 0; i < n; ++i) {
-				const bam_pileup1_t *p = pl + i;
-				if (p->indel == ret->indel1) ++ret->cnt1;
-				else if (p->indel == ret->indel2) ++ret->cnt2;
-				else ++ret->cnt_anti;
-			}
-			{ // write gl[]
-				int tmp, seq_err = 0;
-				double x = 1.0;
-				tmp = max1_i - max2_i;
-				if (tmp < 0) tmp = -tmp;
-				for (j = 0; j < tmp + 1; ++j) x *= INDEL_EXT_DEP;
-				seq_err = mi->q_indel * (1.0 - x) / (1.0 - INDEL_EXT_DEP);
-				ret->gl[0] = ret->gl[1] = 0;
-				for (j = 0; j < n; ++j) {
-					int s1 = pscore[max1_i*n + j], s2 = pscore[max2_i*n + j];
-					//fprintf(stderr, "id=%d, %d, %d, %d, %d, %d\n", j, pl[j].b->core.pos+1, types[max1_i], types[max2_i], s1, s2);
-					if (s1 > s2) ret->gl[0] += s1 - s2 < seq_err? s1 - s2 : seq_err;
-					else ret->gl[1] += s2 - s1 < seq_err? s2 - s1 : seq_err;
-				}
-			}
-			// write cnt_ref and cnt_ambi
-			if (max1_i != 0 && max2_i != 0) {
-				for (j = 0; j < n; ++j) {
-					int diff1 = score[j] - score[max1_i * n + j];
-					int diff2 = score[j] - score[max2_i * n + j];
-					if (diff1 > 0 && diff2 > 0) ++ret->cnt_ref;
-					else if (diff1 == 0 || diff2 == 0) ++ret->cnt_ambi;
-				}
-			}
-		}
-		free(score); free(pscore); free(ref2); free(rs); free(inscns);
-	}
-	{ // call genotype
-		int q[3], qr_indel = (int)(-4.343 * log(mi->r_indel) + 0.5);
-		int min1, min2, min1_i;
-		q[0] = ret->gl[0] + (ret->s[0][0] != '*'? 0 : 0) * qr_indel;
-		q[1] = ret->gl[1] + (ret->s[1][0] != '*'? 0 : 0) * qr_indel;
-		q[2] = n * 3 + (ret->s[0][0] == '*' || ret->s[1][0] == '*'? 1 : 1) * qr_indel;
-		min1 = min2 = 0x7fffffff; min1_i = -1;
-		for (i = 0; i < 3; ++i) {
-			if (q[i] < min1) {
-				min2 = min1; min1 = q[i]; min1_i = i;
-			} else if (q[i] < min2) min2 = q[i];
-		}
-		ret->gt = min1_i;
-		ret->q_cns = min2 - min1;
-		// set q_ref
-		if (ret->gt < 2) ret->q_ref = (ret->s[ret->gt][0] == '*')? 0 : q[1-ret->gt] - q[ret->gt] - qr_indel - 3;
-		else ret->q_ref = (ret->s[0][0] == '*')? q[0] - q[2] : q[1] - q[2];
-		if (ret->q_ref < 0) ret->q_ref = 0;
-	}
-	free(types);
-	return ret;
-}
diff --git a/bam_maqcns.h b/bam_maqcns.h
deleted file mode 100644
index 291ae53..0000000
--- a/bam_maqcns.h
+++ /dev/null
@@ -1,61 +0,0 @@
-#ifndef BAM_MAQCNS_H
-#define BAM_MAQCNS_H
-
-#include "glf.h"
-
-#define BAM_ERRMOD_MAQ2 0
-#define BAM_ERRMOD_MAQ  1
-#define BAM_ERRMOD_SOAP 2
-
-struct __bmc_aux_t;
-
-typedef struct {
-	float het_rate, theta;
-	int n_hap, cap_mapQ, errmod, min_baseQ;
-
-	float eta, q_r;
-	double *fk, *coef;
-	double *lhet;
-	struct __bmc_aux_t *aux;
-} bam_maqcns_t;
-
-typedef struct {
-	int q_indel; // indel sequencing error, phred scaled
-	float r_indel; // indel prior
-	float r_snp; // snp prior
-	// hidden parameters, unchangeable from command line
-	int mm_penalty, indel_err, ambi_thres;
-} bam_maqindel_opt_t;
-
-typedef struct {
-	int indel1, indel2;
-	int cnt1, cnt2, cnt_anti;
-	int cnt_ref, cnt_ambi;
-	char *s[2];
-	//
-	int gt, gl[2];
-	int q_cns, q_ref;
-} bam_maqindel_ret_t;
-
-#ifdef __cplusplus
-extern "C" {
-#endif
-
-	bam_maqcns_t *bam_maqcns_init();
-	void bam_maqcns_prepare(bam_maqcns_t *bm);
-	void bam_maqcns_destroy(bam_maqcns_t *bm);
-	glf1_t *bam_maqcns_glfgen(int n, const bam_pileup1_t *pl, uint8_t ref_base, bam_maqcns_t *bm);
-	uint32_t bam_maqcns_call(int n, const bam_pileup1_t *pl, bam_maqcns_t *bm);
-	// return: cns<<28 | cns2<<24 | mapQ<<16 | cnsQ<<8 | cnsQ2
-	uint32_t glf2cns(const glf1_t *g, int q_r);
-
-	bam_maqindel_opt_t *bam_maqindel_opt_init();
-	bam_maqindel_ret_t *bam_maqindel(int n, int pos, const bam_maqindel_opt_t *mi, const bam_pileup1_t *pl, const char *ref,
-									 int _n_types, int *_types);
-	void bam_maqindel_ret_destroy(bam_maqindel_ret_t*);
-
-#ifdef __cplusplus
-}
-#endif
-
-#endif
diff --git a/bam_plcmd.c b/bam_plcmd.c
index f6381c3..0aa44f5 100644
--- a/bam_plcmd.c
+++ b/bam_plcmd.c
@@ -6,93 +6,8 @@
 #include <errno.h>
 #include "sam.h"
 #include "faidx.h"
-#include "bam_maqcns.h"
-#include "khash.h"
-#include "glf.h"
 #include "kstring.h"
 
-typedef int *indel_list_t;
-KHASH_MAP_INIT_INT64(64, indel_list_t)
-
-#define BAM_PLF_SIMPLE     0x01
-#define BAM_PLF_CNS        0x02
-#define BAM_PLF_INDEL_ONLY 0x04
-#define BAM_PLF_GLF        0x08
-#define BAM_PLF_VAR_ONLY   0x10
-#define BAM_PLF_2ND        0x20
-#define BAM_PLF_RANBASE    0x40
-#define BAM_PLF_1STBASE    0x80
-#define BAM_PLF_ALLBASE    0x100
-#define BAM_PLF_READPOS    0x200
-#define BAM_PLF_NOBAQ      0x400
-
-typedef struct {
-	bam_header_t *h;
-	bam_maqcns_t *c;
-	bam_maqindel_opt_t *ido;
-	faidx_t *fai;
-	khash_t(64) *hash;
-	uint32_t format;
-	int tid, len, last_pos;
-	int mask;
-	int capQ_thres, min_baseQ;
-    int max_depth;  // for indel calling, ignore reads with the depth too high. 0 for unlimited
-	char *ref;
-	glfFile fp_glf; // for glf output only
-} pu_data_t;
-
-char **__bam_get_lines(const char *fn, int *_n);
-void bam_init_header_hash(bam_header_t *header);
-int32_t bam_get_tid(const bam_header_t *header, const char *seq_name);
-
-static khash_t(64) *load_pos(const char *fn, bam_header_t *h)
-{
-	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;
-}
-
 static inline int printw(int c, FILE *fp)
 {
 	char buf[16];
@@ -108,75 +23,6 @@ static inline int printw(int c, FILE *fp)
 	return 0;
 }
 
-// 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);
-	return 0;
-}
-
 static inline void pileup_seq(const bam_pileup1_t *p, int pos, int ref_len, const char *ref)
 {
 	int j;
@@ -212,317 +58,6 @@ static inline void pileup_seq(const bam_pileup1_t *p, int pos, int ref_len, cons
 	if (p->is_tail) putchar('$');
 }
 
-static int pileup_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 i, rb, rms_mapq = -1, *proposed_indels = 0;
-	uint64_t rms_aux;
-	uint32_t cns = 0;
-
-	// if GLF is required, suppress -c completely
-	if (d->format & BAM_PLF_GLF) return glt3_func(tid, pos, n, pu, data);
-	// if d->hash is initialized, only output the sites in the hash table
-	if (d->hash) {
-		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);
-	}
-	// update d->ref if necessary
-	if (d->fai && (int)tid != d->tid) {
-		free(d->ref);
-		d->ref = faidx_fetch_seq(d->fai, d->h->target_name[tid], 0, 0x7fffffff, &d->len);
-		d->tid = tid;
-	}
-	rb = (d->ref && (int)pos < d->len)? d->ref[pos] : 'N';
-	// when the indel-only mode is asked for, return if no reads mapped with indels
-	if (d->format & BAM_PLF_INDEL_ONLY) {
-		for (i = 0; i < n; ++i)
-			if (pu[i].indel != 0) break;
-		if (i == n) return 0;
-	}
-	// call the consensus and indel
-	if (d->format & BAM_PLF_CNS) { // call consensus
-		if (d->format & (BAM_PLF_RANBASE|BAM_PLF_1STBASE)) { // use a random base or the 1st base as the consensus call
-			const bam_pileup1_t *p = (d->format & BAM_PLF_1STBASE)? pu : pu + (int)(drand48() * n);
-			int q = bam1_qual(p->b)[p->qpos];
-			int mapQ = p->b->core.qual < d->c->cap_mapQ? p->b->core.qual : d->c->cap_mapQ;
-			uint32_t b = bam1_seqi(bam1_seq(p->b), p->qpos);
-			cns = b<<28 | 0xf<<24 | mapQ<<16 | q<<8;
-		} else if (d->format & BAM_PLF_ALLBASE) { // collapse all bases
-			uint64_t rmsQ = 0;
-			uint32_t b = 0;
-			for (i = 0; i < n; ++i) {
-				const bam_pileup1_t *p = pu + i;
-				int q = p->b->core.qual < d->c->cap_mapQ? p->b->core.qual : d->c->cap_mapQ;
-				b |= bam1_seqi(bam1_seq(p->b), p->qpos);
-				rmsQ += q * q;
-			}
-			rmsQ = (uint64_t)(sqrt((double)rmsQ / n) + .499);
-			cns = b<<28 | 0xf<<24 | rmsQ<<16 | 60<<8;
-		} else {
-			glf1_t *g = bam_maqcns_glfgen(n, pu, bam_nt16_table[rb], d->c);
-			cns = g->depth == 0? (0xfu<<28 | 0xf<<24) : glf2cns(g, (int)(d->c->q_r + .499));
-			free(g);
-		}
-	}
-    if ((d->format & (BAM_PLF_CNS|BAM_PLF_INDEL_ONLY)) && d->ref && pos < d->len) { // call indels
-        int m = (!d->max_depth || d->max_depth>n) ? n : d->max_depth;
-        if (proposed_indels) // the first element gives the size of the array
-            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);
-	}
-	// when only variant sites are asked for, test if the site is a variant
-	if ((d->format & BAM_PLF_CNS) && (d->format & BAM_PLF_VAR_ONLY)) {
-		if (!(bam_nt16_table[rb] != 15 && cns>>28 != 15 && cns>>28 != bam_nt16_table[rb])) { // not a SNP
-			if (!(r && (r->gt == 2 || strcmp(r->s[r->gt], "*")))) { // not an indel
-				if (r) bam_maqindel_ret_destroy(r);
-				return 0;
-			}
-		}
-	}
-	// print the first 3 columns
-	fputs(d->h->target_name[tid], stdout); putchar('\t');
-	printw(pos+1, stdout); putchar('\t'); putchar(rb); putchar('\t');
-	// print consensus information if required
-	if (d->format & BAM_PLF_CNS) {
-		putchar(bam_nt16_rev_table[cns>>28]); putchar('\t');
-		printw(cns>>8&0xff, stdout); putchar('\t');
-		printw(cns&0xff, stdout); putchar('\t');
-		printw(cns>>16&0xff, stdout); putchar('\t');
-	}
-	// print pileup sequences
-	printw(n, stdout); putchar('\t');
-	for (i = 0; i < n; ++i)
-		pileup_seq(pu + i, pos, d->len, d->ref);
-	// finalize rms_mapq
-	if (d->format & BAM_PLF_CNS) {
-		for (i = rms_aux = 0; i < n; ++i) {
-			const bam_pileup1_t *p = pu + i;
-			int tmp = p->b->core.qual < d->c->cap_mapQ? p->b->core.qual : d->c->cap_mapQ;
-			rms_aux += tmp * tmp;
-		}
-		rms_aux = (uint64_t)(sqrt((double)rms_aux / n) + .499);
-		if (rms_mapq < 0) rms_mapq = rms_aux;
-	}
-	putchar('\t');
-	// print quality
-	for (i = 0; i < n; ++i) {
-		const bam_pileup1_t *p = pu + i;
-		int c = bam1_qual(p->b)[p->qpos] + 33;
-		if (c > 126) c = 126;
-		putchar(c);
-	}
-	if (d->format & BAM_PLF_2ND) { // print 2nd calls and qualities
-		const unsigned char *q;
-		putchar('\t');
-		for (i = 0; i < n; ++i) {
-			const bam_pileup1_t *p = pu + i;
-			q = bam_aux_get(p->b, "E2");
-			putchar(q? q[p->qpos + 1] : 'N');
-		}
-		putchar('\t');
-		for (i = 0; i < n; ++i) {
-			const bam_pileup1_t *p = pu + i;
-			q = bam_aux_get(p->b, "U2");
-			putchar(q? q[p->qpos + 1] : '!');
-		}
-	}
-	// print mapping quality if -s is flagged on the command line
-	if (d->format & BAM_PLF_SIMPLE) {
-		putchar('\t');
-		for (i = 0; i < n; ++i) {
-			int c = pu[i].b->core.qual + 33;
-			if (c > 126) c = 126;
-			putchar(c);
-		}
-	}
-	// print read position
-	if (d->format & BAM_PLF_READPOS) {
-		putchar('\t');
-		for (i = 0; i < n; ++i) {
-			int x = pu[i].qpos;
-			int l = pu[i].b->core.l_qseq;
-			printw(x < l/2? x+1 : -((l-1)-x+1), stdout); putchar(',');
-		}
-	}
-	putchar('\n');
-	// print the indel line if r has been calculated. This only happens if:
-	// a) -c or -i are flagged, AND b) the reference sequence is available
-	if (r) {
-		printf("%s\t%d\t*\t", d->h->target_name[tid], pos + 1);
-		if (r->gt < 2) printf("%s/%s\t", r->s[r->gt], r->s[r->gt]);
-		else printf("%s/%s\t", r->s[0], r->s[1]);
-		printf("%d\t%d\t", r->q_cns, r->q_ref);
-		printf("%d\t%d\t", rms_mapq, n);
-		printf("%s\t%s\t", r->s[0], r->s[1]);
-		//printf("%d\t%d\t", r->gl[0], r->gl[1]);
-		printf("%d\t%d\t%d\t", r->cnt1, r->cnt2, r->cnt_anti);
-		printf("%d\t%d\n", r->cnt_ref, r->cnt_ambi);
-		bam_maqindel_ret_destroy(r);
-	}
-	return 0;
-}
-
-int bam_pileup(int argc, char *argv[])
-{
-	int c, is_SAM = 0;
-	char *fn_list = 0, *fn_fa = 0, *fn_pos = 0;
-	pu_data_t *d = (pu_data_t*)calloc(1, sizeof(pu_data_t));
-    d->max_depth = 1024; d->tid = -1; d->mask = BAM_DEF_MASK; d->min_baseQ = 13;
-	d->c = bam_maqcns_init();
-	d->c->errmod = BAM_ERRMOD_MAQ2; // change the default model
-	d->ido = bam_maqindel_opt_init();
-	while ((c = getopt(argc, argv, "st:f:cT:N:r:l:d:im:gI:G:vM:S2aR:PAQ:C:B")) >= 0) {
-		switch (c) {
-		case 'Q': d->c->min_baseQ = atoi(optarg); break;
-		case 'C': d->capQ_thres = atoi(optarg); break;
-		case 'B': d->format |= BAM_PLF_NOBAQ; break;
-		case 'a': d->c->errmod = BAM_ERRMOD_SOAP; break;
-		case 'A': d->c->errmod = BAM_ERRMOD_MAQ; break;
-		case 's': d->format |= BAM_PLF_SIMPLE; break;
-		case 't': fn_list = strdup(optarg); break;
-		case 'l': fn_pos = strdup(optarg); break;
-		case 'f': fn_fa = strdup(optarg); break;
-		case 'T': d->c->theta = atof(optarg); break;
-		case 'N': d->c->n_hap = atoi(optarg); break;
-		case 'r': d->c->het_rate = atof(optarg); d->ido->r_snp = d->c->het_rate; break;
-		case 'M': d->c->cap_mapQ = atoi(optarg); break;
-		case 'd': d->max_depth = atoi(optarg); break;
-		case 'c': d->format |= BAM_PLF_CNS; break;
-		case 'i': d->format |= BAM_PLF_INDEL_ONLY; break;
-		case 'v': d->format |= BAM_PLF_VAR_ONLY; break;
-		case 'm': d->mask = strtol(optarg, 0, 0); break;
-		case 'g': d->format |= BAM_PLF_GLF; break;
-		case '2': d->format |= BAM_PLF_2ND; break;
-		case 'P': d->format |= BAM_PLF_READPOS; break;
-		case 'I': d->ido->q_indel = atoi(optarg); break;
-		case 'G': d->ido->r_indel = atof(optarg); break;
-		case 'S': is_SAM = 1; break;
-		case 'R':
-			if (strcmp(optarg, "random") == 0) d->format |= BAM_PLF_RANBASE;
-			else if (strcmp(optarg, "first") == 0) d->format |= BAM_PLF_1STBASE;
-			else if (strcmp(optarg, "all") == 0) d->format |= BAM_PLF_ALLBASE;
-			else fprintf(stderr, "[bam_pileup] unrecognized -R\n");
-			break;
-		default: fprintf(stderr, "Unrecognizd option '-%c'.\n", c); return 1;
-		}
-	}
-	if (d->c->errmod != BAM_ERRMOD_MAQ2) d->c->theta += 0.02;
-	if (d->c->theta > 1.0) d->c->theta = 1.0;
-	if (fn_list) is_SAM = 1;
-	if (optind == argc) {
-		fprintf(stderr, "\n");
-		fprintf(stderr, "Usage:  samtools pileup [options] <in.bam>|<in.sam>\n\n");
-		fprintf(stderr, "Option: -s        simple (yet incomplete) pileup format\n");
-		fprintf(stderr, "        -S        the input is in SAM\n");
-		fprintf(stderr, "        -B        disable BAQ computation\n");
-		fprintf(stderr, "        -A        use the original MAQ model for SNP calling (DEPRECATED)\n");
-		fprintf(stderr, "        -2        output the 2nd best call and quality\n");
-		fprintf(stderr, "        -i        only show lines/consensus with indels\n");
-		fprintf(stderr, "        -Q INT    min base quality (possibly capped by BAQ) [%d]\n", d->c->min_baseQ);
-		fprintf(stderr, "        -C INT    coefficient for adjusting mapQ of poor mappings [%d]\n", d->capQ_thres);
-		fprintf(stderr, "        -m INT    filtering reads with bits in INT [0x%x]\n", d->mask);
-		fprintf(stderr, "        -M INT    cap mapping quality at INT [%d]\n", d->c->cap_mapQ);
-        fprintf(stderr, "        -d INT    limit maximum depth for indels [%d]\n", d->max_depth);
-		fprintf(stderr, "        -t FILE   list of reference sequences (force -S)\n");
-		fprintf(stderr, "        -l FILE   list of sites at which pileup is output\n");
-		fprintf(stderr, "        -f FILE   reference sequence in the FASTA format\n\n");
-		fprintf(stderr, "        -c        compute the consensus sequence\n");
-		fprintf(stderr, "        -v        print variants only (for -c)\n");
-		fprintf(stderr, "        -g        output in the GLFv3 format (DEPRECATED)\n");
-		fprintf(stderr, "        -T FLOAT  theta in maq consensus calling model (for -c) [%.4g]\n", d->c->theta);
-		fprintf(stderr, "        -N INT    number of haplotypes in the sample (for -c) [%d]\n", d->c->n_hap);
-		fprintf(stderr, "        -r FLOAT  prior of a difference between two haplotypes (for -c) [%.4g]\n", d->c->het_rate);
-		fprintf(stderr, "        -G FLOAT  prior of an indel between two haplotypes (for -c) [%.4g]\n", d->ido->r_indel);
-		fprintf(stderr, "        -I INT    phred prob. of an indel in sequencing/prep. (for -c) [%d]\n", d->ido->q_indel);
-		fprintf(stderr, "\n");
-		fprintf(stderr, "Warning: Please use the `mpileup' command instead `pileup'. `Pileup' is deprecated!\n\n");
-		free(fn_list); free(fn_fa); free(d);
-		return 1;
-	}
-	if (d->format & (BAM_PLF_RANBASE|BAM_PLF_1STBASE|BAM_PLF_ALLBASE)) d->format |= BAM_PLF_CNS;
-	if (fn_fa) d->fai = fai_load(fn_fa);
-	if (d->format & (BAM_PLF_CNS|BAM_PLF_GLF)) bam_maqcns_prepare(d->c); // consensus calling
-	if (d->format & BAM_PLF_GLF) { // for glf output
-		glf3_header_t *h;
-		h = glf3_header_init();
-		d->fp_glf = bgzf_fdopen(fileno(stdout), "w");
-		glf3_header_write(d->fp_glf, h);
-		glf3_header_destroy(h);
-	}
-	if (d->fai == 0 && (d->format & (BAM_PLF_CNS|BAM_PLF_INDEL_ONLY)))
-		fprintf(stderr, "[bam_pileup] indels will not be called when -f is absent.\n");
-	if (fn_fa && is_SAM && fn_list == 0) fn_list = samfaipath(fn_fa);
-
-	{
-		samfile_t *fp;
-		fp = is_SAM? samopen(argv[optind], "r", fn_list) : samopen(argv[optind], "rb", 0);
-		if (fp == 0 || fp->header == 0) {
-			fprintf(stderr, "[bam_pileup] fail to read the header: non-exisiting file or wrong format.\n");
-			return 1;
-		}
-		d->h = fp->header;
-		if (fn_pos) d->hash = load_pos(fn_pos, d->h);
-		{ // run pileup
-			extern int bam_prob_realn(bam1_t *b, const char *ref);
-			extern int bam_cap_mapQ(bam1_t *b, char *ref, int thres);
-			bam1_t *b;
-			int ret, tid, pos, n_plp;
-			bam_plp_t iter;
-			const bam_pileup1_t *plp;
-			b = bam_init1();
-			iter = bam_plp_init(0, 0);
-			bam_plp_set_mask(iter, d->mask);
-			while ((ret = samread(fp, b)) >= 0) {
-				int skip = 0;
-				if ((int)b->core.tid < 0) break;
-				// update d->ref if necessary
-				if (d->fai && (int)b->core.tid != d->tid) {
-					free(d->ref);
-					d->ref = faidx_fetch_seq(d->fai, d->h->target_name[b->core.tid], 0, 0x7fffffff, &d->len);
-					d->tid = b->core.tid;
-				}
-				if (d->ref && (d->format&BAM_PLF_CNS) && !(d->format&BAM_PLF_NOBAQ)) bam_prob_realn(b, d->ref);
-				if (d->ref && (d->format&BAM_PLF_CNS) && d->capQ_thres > 10) {
-					int q = bam_cap_mapQ(b, d->ref, d->capQ_thres);
-					if (q < 0) skip = 1;
-					else if (b->core.qual > q) b->core.qual = q;
-				} else if (b->core.flag&BAM_FUNMAP) skip = 1;
-				else if ((d->format&BAM_PLF_CNS) && (b->core.flag&1) && !(b->core.flag&2)) skip = 1;
-				if (skip) continue;
-				bam_plp_push(iter, b);
-				while ((plp = bam_plp_next(iter, &tid, &pos, &n_plp)) != 0)
-					pileup_func(tid, pos, n_plp, plp, d);
-			}
-			bam_plp_push(iter, 0);
-			while ((plp = bam_plp_next(iter, &tid, &pos, &n_plp)) != 0)
-				pileup_func(tid, pos, n_plp, plp, d);
-			bam_plp_destroy(iter);
-			bam_destroy1(b);
-		}
-		samclose(fp); // d->h will be destroyed here
-	}
-
-	// free
-	if (d->format & BAM_PLF_GLF) bgzf_close(d->fp_glf);
-	if (fn_pos) { // free the hash table
-		khint_t k;
-		for (k = kh_begin(d->hash); k < kh_end(d->hash); ++k)
-			if (kh_exist(d->hash, k)) free(kh_val(d->hash, k));
-		kh_destroy(64, d->hash);
-	}
-	free(fn_pos); free(fn_list); free(fn_fa);
-	if (d->fai) fai_destroy(d->fai);
-	bam_maqcns_destroy(d->c);
-	free(d->ido); free(d->ref); free(d);
-	return 0;
-}
-
-/***********
- * mpileup *
- ***********/
-
 #include <assert.h>
 #include "bam2bcf.h"
 #include "sample.h"
diff --git a/bam_tview.c b/bam_tview.c
index bf01e15..4eea955 100644
--- a/bam_tview.c
+++ b/bam_tview.c
@@ -19,9 +19,10 @@
 #include <ctype.h>
 #include <assert.h>
 #include <string.h>
+#include <math.h>
 #include "bam.h"
 #include "faidx.h"
-#include "bam_maqcns.h"
+#include "bam2bcf.h"
 
 char bam_aux_getCEi(bam1_t *b, int i);
 char bam_aux_getCSi(bam1_t *b, int i);
@@ -50,7 +51,7 @@ typedef struct {
 	bamFile fp;
 	int curr_tid, left_pos;
 	faidx_t *fai;
-	bam_maqcns_t *bmc;
+	bcf_callaux_t *bca;
 
 	int ccol, last_pos, row_shift, base_for, color_for, is_dot, l_ref, ins, no_skip, show_name;
 	char *ref;
@@ -58,6 +59,7 @@ typedef struct {
 
 int tv_pl_func(uint32_t tid, uint32_t pos, int n, const bam_pileup1_t *pl, void *data)
 {
+	extern unsigned char bam_nt16_table[256];
 	tview_t *tv = (tview_t*)data;
 	int i, j, c, rb, attr, max_ins = 0;
 	uint32_t call = 0;
@@ -70,11 +72,26 @@ int tv_pl_func(uint32_t tid, uint32_t pos, int n, const bam_pileup1_t *pl, void
 		mvaddch(1, tv->ccol++, c);
 	}
 	if (pos%10 == 0 && tv->mcol - tv->ccol >= 10) mvprintw(0, tv->ccol, "%-d", pos+1);
-	// print consensus
-	call = bam_maqcns_call(n, pl, tv->bmc);
+	{ // call consensus
+		bcf_callret1_t bcr;
+		int qsum[4], a1, a2, tmp;
+		double p[3], prior = 30;
+		bcf_call_glfgen(n, pl, bam_nt16_table[rb], tv->bca, &bcr);
+		for (i = 0; i < 4; ++i) qsum[i] = bcr.qsum[i]<<2 | i;
+		for (i = 1; i < 4; ++i) // insertion sort
+			for (j = i; j > 0 && qsum[j] > qsum[j-1]; --j)
+				tmp = qsum[j], qsum[j] = qsum[j-1], qsum[j-1] = tmp;
+		a1 = qsum[0]&3; a2 = qsum[1]&3;
+		p[0] = bcr.p[a1*5+a1]; p[1] = bcr.p[a1*5+a2] + prior; p[2] = bcr.p[a2*5+a2];
+		if ("ACGT"[a1] != toupper(rb)) p[0] += prior + 3;
+		if ("ACGT"[a2] != toupper(rb)) p[2] += prior + 3;
+		if (p[0] < p[1] && p[0] < p[2]) call = (1<<a1)<<16 | (int)((p[1]<p[2]?p[1]:p[2]) - p[0] + .499);
+		else if (p[2] < p[1] && p[2] < p[0]) call = (1<<a2)<<16 | (int)((p[0]<p[1]?p[0]:p[1]) - p[2] + .499);
+		else call = (1<<a1|1<<a2)<<16 | (int)((p[0]<p[2]?p[0]:p[2]) - p[1] + .499);
+	}
 	attr = A_UNDERLINE;
-	c = ",ACMGRSVTWYHKDBN"[call>>28&0xf];
-	i = (call>>8&0xff)/10+1;
+	c = ",ACMGRSVTWYHKDBN"[call>>16&0xf];
+	i = (call&0xffff)/10+1;
 	if (i > 4) i = 4;
 	attr |= COLOR_PAIR(i);
 	if (c == toupper(rb)) c = '.';
@@ -191,9 +208,8 @@ tview_t *tv_init(const char *fn, const char *fn_fa)
 	if (tv->idx == 0) exit(1);
 	tv->lplbuf = bam_lplbuf_init(tv_pl_func, tv);
 	if (fn_fa) tv->fai = fai_load(fn_fa);
-	tv->bmc = bam_maqcns_init();
+	tv->bca = bcf_call_init(0.83, 13);
 	tv->ins = 1;
-	bam_maqcns_prepare(tv->bmc);
 
 	initscr();
 	keypad(stdscr, TRUE);
@@ -224,7 +240,7 @@ void tv_destroy(tview_t *tv)
 	endwin();
 
 	bam_lplbuf_destroy(tv->lplbuf);
-	bam_maqcns_destroy(tv->bmc);
+	bcf_call_destroy(tv->bca);
 	bam_index_destroy(tv->idx);
 	if (tv->fai) fai_destroy(tv->fai);
 	free(tv->ref);
diff --git a/bamtk.c b/bamtk.c
index 0941f6f..79b0dbc 100644
--- a/bamtk.c
+++ b/bamtk.c
@@ -9,7 +9,6 @@
 #endif
 
 int bam_taf2baf(int argc, char *argv[]);
-int bam_pileup(int argc, char *argv[]);
 int bam_mpileup(int argc, char *argv[]);
 int bam_merge(int argc, char *argv[]);
 int bam_index(int argc, char *argv[]);
@@ -29,7 +28,6 @@ int main_cat(int argc, char *argv[]);
 int main_depth(int argc, char *argv[]);
 
 int faidx_main(int argc, char *argv[]);
-int glf3_view_main(int argc, char *argv[]);
 
 static int usage()
 {
@@ -39,7 +37,6 @@ static int usage()
 	fprintf(stderr, "Usage:   samtools <command> [options]\n\n");
 	fprintf(stderr, "Command: view        SAM<->BAM conversion\n");
 	fprintf(stderr, "         sort        sort alignment file\n");
-	fprintf(stderr, "         pileup      generate pileup output\n");
 	fprintf(stderr, "         mpileup     multi-way pileup\n");
 	fprintf(stderr, "         depth       compute the depth\n");
 	fprintf(stderr, "         faidx       index/extract FASTA\n");
@@ -49,7 +46,6 @@ static int usage()
 	fprintf(stderr, "         index       index alignment\n");
 	fprintf(stderr, "         idxstats    BAM index stats (r595 or later)\n");
 	fprintf(stderr, "         fixmate     fix mate information\n");
-	fprintf(stderr, "         glfview     print GLFv3 file\n");
 	fprintf(stderr, "         flagstat    simple stats\n");
 	fprintf(stderr, "         calmd       recalculate MD/NM tags and '=' bases\n");
 	fprintf(stderr, "         merge       merge sorted alignments\n");
@@ -80,7 +76,6 @@ int main(int argc, char *argv[])
 	if (argc < 2) return usage();
 	if (strcmp(argv[1], "view") == 0) return main_samview(argc-1, argv+1);
 	else if (strcmp(argv[1], "import") == 0) return main_import(argc-1, argv+1);
-	else if (strcmp(argv[1], "pileup") == 0) return bam_pileup(argc-1, argv+1);
 	else if (strcmp(argv[1], "mpileup") == 0) return bam_mpileup(argc-1, argv+1);
 	else if (strcmp(argv[1], "merge") == 0) return bam_merge(argc-1, argv+1);
 	else if (strcmp(argv[1], "sort") == 0) return bam_sort(argc-1, argv+1);
@@ -89,7 +84,6 @@ int main(int argc, char *argv[])
 	else if (strcmp(argv[1], "faidx") == 0) return faidx_main(argc-1, argv+1);
 	else if (strcmp(argv[1], "fixmate") == 0) return bam_mating(argc-1, argv+1);
 	else if (strcmp(argv[1], "rmdup") == 0) return bam_rmdup(argc-1, argv+1);
-	else if (strcmp(argv[1], "glfview") == 0) return glf3_view_main(argc-1, argv+1);
 	else if (strcmp(argv[1], "flagstat") == 0) return bam_flagstat(argc-1, argv+1);
 	else if (strcmp(argv[1], "calmd") == 0) return bam_fillmd(argc-1, argv+1);
 	else if (strcmp(argv[1], "fillmd") == 0) return bam_fillmd(argc-1, argv+1);
diff --git a/glf.c b/glf.c
deleted file mode 100644
index 8d5346a..0000000
--- a/glf.c
+++ /dev/null
@@ -1,236 +0,0 @@
-#include <string.h>
-#include <stdlib.h>
-#include "glf.h"
-
-#ifdef _NO_BGZF
-// then alias bgzf_*() functions
-#endif
-
-static int glf3_is_BE = 0;
-
-static inline uint32_t bam_swap_endian_4(uint32_t v)
-{
-	v = ((v & 0x0000FFFFU) << 16) | (v >> 16);
-	return ((v & 0x00FF00FFU) << 8) | ((v & 0xFF00FF00U) >> 8);
-}
-
-static inline uint16_t bam_swap_endian_2(uint16_t v)
-{
-	return (uint16_t)(((v & 0x00FF00FFU) << 8) | ((v & 0xFF00FF00U) >> 8));
-}
-
-static inline int bam_is_big_endian()
-{
-	long one= 1;
-	return !(*((char *)(&one)));
-}
-
-glf3_header_t *glf3_header_init()
-{
-	glf3_is_BE = bam_is_big_endian();
-	return (glf3_header_t*)calloc(1, sizeof(glf3_header_t));
-}
-
-glf3_header_t *glf3_header_read(glfFile fp)
-{
-	glf3_header_t *h;
-	char magic[4];
-	h = glf3_header_init();
-	bgzf_read(fp, magic, 4);
-	if (strncmp(magic, "GLF\3", 4)) {
-		fprintf(stderr, "[glf3_header_read] invalid magic.\n");
-		glf3_header_destroy(h);
-		return 0;
-	}
-	bgzf_read(fp, &h->l_text, 4);
-	if (glf3_is_BE) h->l_text = bam_swap_endian_4(h->l_text);
-	if (h->l_text) {
-		h->text = (uint8_t*)calloc(h->l_text + 1, 1);
-		bgzf_read(fp, h->text, h->l_text);
-	}
-	return h;
-}
-
-void glf3_header_write(glfFile fp, const glf3_header_t *h)
-{
-	int32_t x;
-	bgzf_write(fp, "GLF\3", 4);
-	x = glf3_is_BE? bam_swap_endian_4(h->l_text) : h->l_text;
-	bgzf_write(fp, &x, 4);
-	if (h->l_text) bgzf_write(fp, h->text, h->l_text);
-}
-
-void glf3_header_destroy(glf3_header_t *h)
-{
-	free(h->text);
-	free(h);
-}
-
-char *glf3_ref_read(glfFile fp, int *len)
-{
-	int32_t n, x;
-	char *str;
-	*len = 0;
-	if (bgzf_read(fp, &n, 4) != 4) return 0;
-	if (glf3_is_BE) n = bam_swap_endian_4(n);
-	if (n < 0) {
-		fprintf(stderr, "[glf3_ref_read] invalid reference name length: %d.\n", n);
-		return 0;
-	}
-	str = (char*)calloc(n + 1, 1); // not necesarily n+1 in fact
-	x = bgzf_read(fp, str, n);
-	x += bgzf_read(fp, len, 4);
-	if (x != n + 4) {
-		free(str); *len = -1; return 0; // truncated
-	}
-	if (glf3_is_BE) *len = bam_swap_endian_4(*len);
-	return str;
-}
-
-void glf3_ref_write(glfFile fp, const char *str, int len)
-{
-	int32_t m, n = strlen(str) + 1;
-	m = glf3_is_BE? bam_swap_endian_4(n) : n;
-	bgzf_write(fp, &m, 4);
-	bgzf_write(fp, str, n);
-	if (glf3_is_BE) len = bam_swap_endian_4(len);
-	bgzf_write(fp, &len, 4);
-}
-
-void glf3_view1(const char *ref_name, const glf3_t *g3, int pos)
-{
-	int j;
-	if (g3->rtype == GLF3_RTYPE_END) return;
-	printf("%s\t%d\t%c\t%d\t%d\t%d", ref_name, pos + 1,
-		   g3->rtype == GLF3_RTYPE_INDEL? '*' : "XACMGRSVTWYHKDBN"[g3->ref_base],
-		   g3->depth, g3->rms_mapQ, g3->min_lk);
-	if (g3->rtype == GLF3_RTYPE_SUB)
-		for (j = 0; j != 10; ++j) printf("\t%d", g3->lk[j]);
-	else {
-		printf("\t%d\t%d\t%d\t%d\t%d\t%s\t%s\t", g3->lk[0], g3->lk[1], g3->lk[2], g3->indel_len[0], g3->indel_len[1],
-			   g3->indel_len[0]? g3->indel_seq[0] : "*", g3->indel_len[1]? g3->indel_seq[1] : "*");
-	}
-	printf("\n");
-}
-
-int glf3_write1(glfFile fp, const glf3_t *g3)
-{
-	int r;
-	uint8_t c;
-	uint32_t y[2];
-	c = g3->rtype<<4 | g3->ref_base;
-	r = bgzf_write(fp, &c, 1);
-	if (g3->rtype == GLF3_RTYPE_END) return r;
-	y[0] = g3->offset;
-	y[1] = g3->min_lk<<24 | g3->depth;
-	if (glf3_is_BE) {
-		y[0] = bam_swap_endian_4(y[0]);
-		y[1] = bam_swap_endian_4(y[1]);
-	}
-	r += bgzf_write(fp, y, 8);
-	r += bgzf_write(fp, &g3->rms_mapQ, 1);
-	if (g3->rtype == GLF3_RTYPE_SUB) r += bgzf_write(fp, g3->lk, 10);
-	else {
-		int16_t x[2];
-		r += bgzf_write(fp, g3->lk, 3);
-		x[0] = glf3_is_BE? bam_swap_endian_2(g3->indel_len[0]) : g3->indel_len[0];
-		x[1] = glf3_is_BE? bam_swap_endian_2(g3->indel_len[1]) : g3->indel_len[1];
-		r += bgzf_write(fp, x, 4);
-		if (g3->indel_len[0]) r += bgzf_write(fp, g3->indel_seq[0], abs(g3->indel_len[0]));
-		if (g3->indel_len[1]) r += bgzf_write(fp, g3->indel_seq[1], abs(g3->indel_len[1]));
-	}
-	return r;
-}
-
-#ifndef kv_roundup32
-#define kv_roundup32(x) (--(x), (x)|=(x)>>1, (x)|=(x)>>2, (x)|=(x)>>4, (x)|=(x)>>8, (x)|=(x)>>16, ++(x))
-#endif
-
-int glf3_read1(glfFile fp, glf3_t *g3)
-{
-	int r;
-	uint8_t c;
-	uint32_t y[2];
-	r = bgzf_read(fp, &c, 1);
-	if (r == 0) return 0;
-	g3->ref_base = c & 0xf;
-	g3->rtype = c>>4;
-	if (g3->rtype == GLF3_RTYPE_END) return r;
-	r += bgzf_read(fp, y, 8);
-	if (glf3_is_BE) {
-		y[0] = bam_swap_endian_4(y[0]);
-		y[1] = bam_swap_endian_4(y[1]);
-	}
-	g3->offset = y[0];
-	g3->min_lk = y[1]>>24;
-	g3->depth = y[1]<<8>>8;
-	r += bgzf_read(fp, &g3->rms_mapQ, 1);
-	if (g3->rtype == GLF3_RTYPE_SUB) r += bgzf_read(fp, g3->lk, 10);
-	else {
-		int16_t x[2], max;
-		r += bgzf_read(fp, g3->lk, 3);
-		r += bgzf_read(fp, x, 4);
-		if (glf3_is_BE) {
-			x[0] = bam_swap_endian_2(x[0]);
-			x[1] = bam_swap_endian_2(x[1]);
-		}
-		g3->indel_len[0] = x[0];
-		g3->indel_len[1] = x[1];
-		x[0] = abs(x[0]); x[1] = abs(x[1]);
-		max = (x[0] > x[1]? x[0] : x[1]) + 1;
-		if (g3->max_len < max) {
-			g3->max_len = max;
-			kv_roundup32(g3->max_len);
-			g3->indel_seq[0] = (char*)realloc(g3->indel_seq[0], g3->max_len);
-			g3->indel_seq[1] = (char*)realloc(g3->indel_seq[1], g3->max_len);
-		}
-		r += bgzf_read(fp, g3->indel_seq[0], x[0]);
-		r += bgzf_read(fp, g3->indel_seq[1], x[1]);
-		g3->indel_seq[0][x[0]] = g3->indel_seq[1][x[1]] = 0;
-	}
-	return r;
-}
-
-void glf3_view(glfFile fp)
-{
-	glf3_header_t *h;
-	char *name;
-	glf3_t *g3;
-	int len;
-	h = glf3_header_read(fp);
-	g3 = glf3_init1();
-	while ((name = glf3_ref_read(fp, &len)) != 0) {
-		int pos = 0;
-		while (glf3_read1(fp, g3) && g3->rtype != GLF3_RTYPE_END) {
-			pos += g3->offset;
-			glf3_view1(name, g3, pos);
-		}
-		free(name);
-	}
-	glf3_header_destroy(h);
-	glf3_destroy1(g3);
-}
-
-int glf3_view_main(int argc, char *argv[])
-{
-	glfFile fp;
-	if (argc == 1) {
-		fprintf(stderr, "Usage: glfview <in.glf>\n");
-		return 1;
-	}
-	fp = (strcmp(argv[1], "-") == 0)? bgzf_fdopen(fileno(stdin), "r") : bgzf_open(argv[1], "r");
-	if (fp == 0) {
-		fprintf(stderr, "Fail to open file '%s'\n", argv[1]);
-		return 1;
-	}
-	glf3_view(fp);
-	bgzf_close(fp);
-	return 0;
-}
-
-#ifdef GLFVIEW_MAIN
-int main(int argc, char *argv[])
-{
-	return glf3_view_main(argc, argv);
-}
-#endif
diff --git a/glf.h b/glf.h
deleted file mode 100644
index 12e5400..0000000
--- a/glf.h
+++ /dev/null
@@ -1,56 +0,0 @@
-#ifndef GLF_H_
-#define GLF_H_
-
-typedef struct {
-	unsigned char ref_base:4, dummy:4; /** "XACMGRSVTWYHKDBN"[ref_base] gives the reference base */
-	unsigned char max_mapQ; /** maximum mapping quality */
-	unsigned char lk[10];   /** log likelihood ratio, capped at 255 */
-	unsigned min_lk:8, depth:24; /** minimum lk capped at 255, and the number of mapped reads */
-} glf1_t;
-
-#include <stdint.h>
-#include "bgzf.h"
-typedef BGZF *glfFile;
-
-#define GLF3_RTYPE_END   0
-#define GLF3_RTYPE_SUB   1
-#define GLF3_RTYPE_INDEL 2
-
-typedef struct {
-	uint8_t ref_base:4, rtype:4; /** "XACMGRSVTWYHKDBN"[ref_base] gives the reference base */
-	uint8_t rms_mapQ; /** RMS mapping quality */
-	uint8_t lk[10];   /** log likelihood ratio, capped at 255 */
-	uint32_t min_lk:8, depth:24; /** minimum lk capped at 255, and the number of mapped reads */
-	int32_t offset; /** the first base in a chromosome has offset zero. */
-	// for indel (lkHom1, lkHom2 and lkHet are the first three elements in lk[10])
-	int16_t indel_len[2];
-	int32_t max_len; // maximum indel len; will be modified by glf3_read1()
-	char *indel_seq[2];
-} glf3_t;
-
-typedef struct {
-	int32_t l_text;
-	uint8_t *text;
-} glf3_header_t;
-
-#ifdef __cplusplus
-extern "C" {
-#endif
-
-#define glf3_init1() ((glf3_t*)calloc(1, sizeof(glf3_t)))
-#define glf3_destroy1(g3) do { free((g3)->indel_seq[0]); free((g3)->indel_seq[1]); free(g3); } while (0)
-
-	glf3_header_t *glf3_header_init();
-	glf3_header_t *glf3_header_read(glfFile fp);
-	void glf3_header_write(glfFile fp, const glf3_header_t *h);
-	void glf3_header_destroy(glf3_header_t *h);
-	char *glf3_ref_read(glfFile fp, int *len);
-	void glf3_ref_write(glfFile fp, const char *name, int len);
-	int glf3_write1(glfFile fp, const glf3_t *g3);
-	int glf3_read1(glfFile fp, glf3_t *g3);
-
-#ifdef __cplusplus
-}
-#endif
-
-#endif
-- 
2.39.2