annotate PsiCLASS-1.0.2/samtools-0.1.19/bcftools/bcf2qcall.c @ 0:903fc43d6227 draft default tip

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