From 7c51db2e81d77efe6526803fa8d88fead5c82131 Mon Sep 17 00:00:00 2001 From: Heng Li Date: Wed, 15 Sep 2010 04:18:31 +0000 Subject: [PATCH] * convert BCF to QCALL input --- bcftools/Makefile | 2 +- bcftools/bcf2qcall.c | 91 ++++++++++++++++++++++++++++++++++++++++++++ bcftools/call1.c | 16 ++++++-- 3 files changed, 104 insertions(+), 5 deletions(-) create mode 100644 bcftools/bcf2qcall.c diff --git a/bcftools/Makefile b/bcftools/Makefile index 5a3dc76..9592025 100644 --- a/bcftools/Makefile +++ b/bcftools/Makefile @@ -1,7 +1,7 @@ CC= gcc CFLAGS= -g -Wall -O2 #-m64 #-arch ppc DFLAGS= -D_FILE_OFFSET_BITS=64 -D_USE_KNETFILE -LOBJS= bcf.o vcf.o bcfutils.o prob1.o kfunc.o index.o fet.o +LOBJS= bcf.o vcf.o bcfutils.o prob1.o kfunc.o index.o fet.o bcf2qcall.o OMISC= .. AOBJS= call1.o main.o $(OMISC)/kstring.o $(OMISC)/bgzf.o $(OMISC)/knetfile.o PROG= bcftools diff --git a/bcftools/bcf2qcall.c b/bcftools/bcf2qcall.c new file mode 100644 index 0000000..8634c9e --- /dev/null +++ b/bcftools/bcf2qcall.c @@ -0,0 +1,91 @@ +#include +#include +#include +#include +#include "bcf.h" + +static int8_t nt4_table[256] = { + 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, + 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, + 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4 /*'-'*/, 4, 4, + 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, + 4, 0, 4, 1, 4, 4, 4, 2, 4, 4, 4, 4, 4, 4, 4, 4, + 4, 4, 4, 4, 3, 4, 4, 4, -1, 4, 4, 4, 4, 4, 4, 4, + 4, 0, 4, 1, 4, 4, 4, 2, 4, 4, 4, 4, 4, 4, 4, 4, + 4, 4, 4, 4, 3, 4, 4, 4, -1, 4, 4, 4, 4, 4, 4, 4, + 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, + 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, + 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, + 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, + 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, + 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, + 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, + 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4 +}; + +static int read_I16(bcf1_t *b, int anno[16]) +{ + char *p; + int i; + if ((p = strstr(b->info, "I16=")) == 0) return -1; + p += 4; + for (i = 0; i < 16; ++i) { + anno[i] = strtol(p, &p, 10); + if (anno[i] == 0 && (errno == EINVAL || errno == ERANGE)) return -2; + ++p; + } + return 0; +} + +int bcf_2qcall(bcf_hdr_t *h, bcf1_t *b) +{ + int a[4], k, g[10], l, map[4], k1, j, i, i0, anno[16], dp, mq, d_rest; + char *s; + if (b->ref[1] != 0 || b->n_alleles > 4) return -1; // ref is not a single base + for (i = 0; i < b->n_gi; ++i) + if (b->gi[i].fmt == bcf_str2int("PL", 2)) break; + if (i == b->n_gi) return -1; // no PL + if (read_I16(b, anno) != 0) return -1; // no I16; FIXME: can be improved + d_rest = dp = anno[0] + anno[1] + anno[2] + anno[3]; + if (dp == 0) return -1; // depth is zero + mq = (int)(sqrt((double)(anno[9] + anno[11]) / dp) + .499); + i0 = i; + a[0] = nt4_table[(int)b->ref[0]]; + if (a[0] > 3) return -1; // ref is not A/C/G/T + a[1] = a[2] = a[3] = -2; // -1 has a special meaning + if (b->alt[0] == 0) return -1; // no alternate allele + map[0] = map[1] = map[2] = map[3] = -2; + map[a[0]] = 0; + for (k = 0, s = b->alt, k1 = -1; k < 3 && *s; ++k, s += 2) { + if (s[1] != ',' && s[1] != 0) return -1; // ALT is not single base + a[k+1] = nt4_table[(int)*s]; + if (a[k+1] >= 0) map[a[k+1]] = k+1; + else k1 = k+1; + if (s[1] == 0) break; + } + for (k = 0; k < 4; ++k) + if (map[k] < 0) map[k] = k1; + for (i = 0; i < h->n_smpl; ++i) { + int d; + uint8_t *p = b->gi[i0].data + i * b->gi[i0].len; + for (j = 0; j < b->gi[i0].len; ++j) + if (p[j]) break; + d = (int)((double)d_rest / (h->n_smpl - i) + .499); + if (d == 0) d = 1; + if (j == b->gi[i0].len) d = 0; + d_rest -= d; + for (k = j = 0; k < 4; ++k) { + for (l = k; l < 4; ++l) { + int t, x = map[k], y = map[l]; + if (x > y) t = x, x = y, y = t; + g[j++] = p[x * b->n_alleles - x * (x-1) / 2 + (y - x)]; + } + } + printf("%s\t%d\t%c", h->ns[b->tid], b->pos+1, *b->ref); + printf("\t%d\t%d\t0", d, mq); + for (j = 0; j < 10; ++j) + printf("\t%d", g[j]); + printf("\t%s\n", h->sns[i]); + } + return 0; +} diff --git a/bcftools/call1.c b/bcftools/call1.c index c0c24a6..08d759d 100644 --- a/bcftools/call1.c +++ b/bcftools/call1.c @@ -23,6 +23,7 @@ KSTREAM_INIT(gzFile, gzread, 16384) #define VC_HWE 128 #define VC_KEEPALT 256 #define VC_ACGT_ONLY 512 +#define VC_QCALL 1024 typedef struct { int flag, prior_type, n1; @@ -182,6 +183,7 @@ static int update_bcf1(int n_smpl, bcf1_t *b, const bcf_p1aux_t *pa, const bcf_p int bcfview(int argc, char *argv[]) { + extern int bcf_2qcall(bcf_hdr_t *h, bcf1_t *b); bcf_t *bp, *bout = 0; bcf1_t *b; int c; @@ -196,7 +198,7 @@ int bcfview(int argc, char *argv[]) tid = begin = end = -1; memset(&vc, 0, sizeof(viewconf_t)); vc.prior_type = vc.n1 = -1; vc.theta = 1e-3; vc.pref = 0.5; - while ((c = getopt(argc, argv, "N1:l:cHAGvLbSuP:t:p:")) >= 0) { + while ((c = getopt(argc, argv, "N1:l:cHAGvLbSuP:t:p:Q")) >= 0) { switch (c) { case '1': vc.n1 = atoi(optarg); break; case 'l': vc.fn_list = strdup(optarg); break; @@ -212,6 +214,7 @@ int bcfview(int argc, char *argv[]) case 'H': vc.flag |= VC_HWE; break; case 't': vc.theta = atof(optarg); break; case 'p': vc.pref = atof(optarg); break; + case 'Q': vc.flag |= VC_QCALL; break; case 'P': if (strcmp(optarg, "full") == 0) vc.prior_type = MC_PTYPE_FULL; else if (strcmp(optarg, "cond2") == 0) vc.prior_type = MC_PTYPE_COND2; @@ -233,6 +236,7 @@ int bcfview(int argc, char *argv[]) fprintf(stderr, " -H perform Hardy-Weinberg test (slower)\n"); fprintf(stderr, " -v output potential variant sites only\n"); fprintf(stderr, " -N skip sites where REF is not A/C/G/T\n"); + fprintf(stderr, " -Q output the QCALL likelihood format\n"); fprintf(stderr, " -1 INT number of group-1 samples [0]\n"); fprintf(stderr, " -l FILE list of sites to output [all sites]\n"); fprintf(stderr, " -t FLOAT scaled mutation rate [%.4lg]\n", vc.theta); @@ -251,7 +255,7 @@ int bcfview(int argc, char *argv[]) bp = vcf_open(argv[optind], moder); h = vcf_hdr_read(bp); bout = vcf_open("-", modew); - vcf_hdr_write(bout, h); + if (!(vc.flag & VC_QCALL)) vcf_hdr_write(bout, h); if (vc.flag & VC_CALL) { p1 = bcf_p1_init(h->n_smpl); if (vc.prior_file) { @@ -300,7 +304,11 @@ int bcfview(int argc, char *argv[]) if (!(l > begin && end > b->pos)) continue; } ++n_processed; - if (vc.flag & VC_CALL) { + if (vc.flag & VC_QCALL) { // output QCALL format; STOP here + bcf_2qcall(h, b); + continue; + } + if (vc.flag & VC_CALL) { // call variants bcf_p1rst_t pr; bcf_gl2pl(h->n_smpl, b); bcf_p1_cal(b, p1, &pr); // pr.g[3] is not calculated here @@ -312,7 +320,7 @@ int bcfview(int argc, char *argv[]) if (pr.p_ref >= vc.pref && (vc.flag & VC_VARONLY)) continue; update_bcf1(h->n_smpl, b, p1, &pr, vc.pref, vc.flag); } - if (vc.flag & VC_NO_GENO) { + if (vc.flag & VC_NO_GENO) { // do not output GENO fields b->n_gi = 0; b->fmt[0] = '\0'; } -- 2.39.5