--- /dev/null
+#include <errno.h>
+#include <math.h>
+#include <string.h>
+#include <stdlib.h>
+#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;
+}
#define VC_HWE 128
#define VC_KEEPALT 256
#define VC_ACGT_ONLY 512
+#define VC_QCALL 1024
typedef struct {
int flag, prior_type, n1;
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;
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;
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;
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);
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) {
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
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';
}