samtools/bcftools/bcf2qcall.c.pysam.c

   1 #include "pysam.h"
   2
   3 #include <errno.h>
   4 #include <math.h>
   5 #include <string.h>
   6 #include <stdlib.h>
   7 #include "bcf.h"
   8
   9 static int8_t nt4_table[256] = {
  10         4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,
  11         4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,
  12         4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,  4, 4 /*'-'*/, 4, 4,
  13         4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,
  14         4, 0, 4, 1,  4, 4, 4, 2,  4, 4, 4, 4,  4, 4, 4, 4,
  15         4, 4, 4, 4,  3, 4, 4, 4, -1, 4, 4, 4,  4, 4, 4, 4,
  16         4, 0, 4, 1,  4, 4, 4, 2,  4, 4, 4, 4,  4, 4, 4, 4,
  17         4, 4, 4, 4,  3, 4, 4, 4, -1, 4, 4, 4,  4, 4, 4, 4,
  18         4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,
  19         4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,
  20         4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,
  21         4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,
  22         4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,
  23         4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,
  24         4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,
  25         4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4
  26 };
  27
  28 static int read_I16(bcf1_t *b, int anno[16])
  29 {
  30         char *p;
  31         int i;
  32         if ((p = strstr(b->info, "I16=")) == 0) return -1;
  33         p += 4;
  34         for (i = 0; i < 16; ++i) {
  35                 anno[i] = strtol(p, &p, 10);
  36                 if (anno[i] == 0 && (errno == EINVAL || errno == ERANGE)) return -2;
  37                 ++p;
  38         }
  39         return 0;
  40 }
  41
  42 int bcf_2qcall(bcf_hdr_t *h, bcf1_t *b)
  43 {
  44         int a[4], k, g[10], l, map[4], k1, j, i, i0, anno[16], dp, mq, d_rest;
  45         char *s;
  46         if (b->ref[1] != 0 || b->n_alleles > 4) return -1; // ref is not a single base
  47         for (i = 0; i < b->n_gi; ++i)
  48                 if (b->gi[i].fmt == bcf_str2int("PL", 2)) break;
  49         if (i == b->n_gi) return -1; // no PL
  50         if (read_I16(b, anno) != 0) return -1; // no I16; FIXME: can be improved
  51         d_rest = dp = anno[0] + anno[1] + anno[2] + anno[3];
  52         if (dp == 0) return -1; // depth is zero
  53         mq = (int)(sqrt((double)(anno[9] + anno[11]) / dp) + .499);
  54         i0 = i;
  55         a[0] = nt4_table[(int)b->ref[0]];
  56         if (a[0] > 3) return -1; // ref is not A/C/G/T
  57         a[1] = a[2] = a[3] = -2; // -1 has a special meaning
  58         if (b->alt[0] == 0) return -1; // no alternate allele
  59         map[0] = map[1] = map[2] = map[3] = -2;
  60         map[a[0]] = 0;
  61         for (k = 0, s = b->alt, k1 = -1; k < 3 && *s; ++k, s += 2) {
  62                 if (s[1] != ',' && s[1] != 0) return -1; // ALT is not single base
  63                 a[k+1] = nt4_table[(int)*s];
  64                 if (a[k+1] >= 0) map[a[k+1]] = k+1;
  65                 else k1 = k+1;
  66                 if (s[1] == 0) break;
  67         }
  68         for (k = 0; k < 4; ++k)
  69                 if (map[k] < 0) map[k] = k1;
  70         for (i = 0; i < h->n_smpl; ++i) {
  71                 int d;
  72                 uint8_t *p = b->gi[i0].data + i * b->gi[i0].len;
  73                 for (j = 0; j < b->gi[i0].len; ++j)
  74                         if (p[j]) break;
  75                 d = (int)((double)d_rest / (h->n_smpl - i) + .499);
  76                 if (d == 0) d = 1;
  77                 if (j == b->gi[i0].len) d = 0;
  78                 d_rest -= d;
  79                 for (k = j = 0; k < 4; ++k) {
  80                         for (l = k; l < 4; ++l) {
  81                                 int t, x = map[k], y = map[l];
  82                                 if (x > y) t = x, x = y, y = t; // swap
  83                                 g[j++] = p[y * (y+1) / 2 + x];
  84                         }
  85                 }
  86                 printf("%s\t%d\t%c", h->ns[b->tid], b->pos+1, *b->ref);
  87                 printf("\t%d\t%d\t0", d, mq);
  88                 for (j = 0; j < 10; ++j)
  89                         printf("\t%d", g[j]);
  90                 printf("\t%s\n", h->sns[i]);
  91         }
  92         return 0;
  93 }