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
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
#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
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;
+++ /dev/null
-#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;
-}
+++ /dev/null
-#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
#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];
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;
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"
#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);
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;
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;
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 = '.';
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);
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);
#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[]);
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()
{
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");
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");
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);
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);
+++ /dev/null
-#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
+++ /dev/null
-#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
projects / samtools.git / commitdiff