|
0
|
1 #include <stdlib.h>
|
|
|
2 #include <stdint.h>
|
|
|
3 #include "bcf.h"
|
|
|
4
|
|
|
5 #define MAX_GENO 359
|
|
|
6
|
|
|
7 int8_t seq_bitcnt[] = { 4, 1, 1, 2, 1, 2, 2, 3, 1, 2, 2, 3, 2, 3, 3, 4 };
|
|
|
8 char *seq_nt16rev = "XACMGRSVTWYHKDBN";
|
|
|
9
|
|
|
10 uint32_t *bcf_trio_prep(int is_x, int is_son)
|
|
|
11 {
|
|
|
12 int i, j, k, n, map[10];
|
|
|
13 uint32_t *ret;
|
|
|
14 ret = calloc(MAX_GENO, 4);
|
|
|
15 for (i = 0, k = 0; i < 4; ++i)
|
|
|
16 for (j = i; j < 4; ++j)
|
|
|
17 map[k++] = 1<<i|1<<j;
|
|
|
18 for (i = 0, n = 1; i < 10; ++i) { // father
|
|
|
19 if (is_x && seq_bitcnt[map[i]] != 1) continue;
|
|
|
20 if (is_x && is_son) {
|
|
|
21 for (j = 0; j < 10; ++j) // mother
|
|
|
22 for (k = 0; k < 10; ++k) // child
|
|
|
23 if (seq_bitcnt[map[k]] == 1 && (map[j]&map[k]))
|
|
|
24 ret[n++] = j<<16 | i<<8 | k;
|
|
|
25 } else {
|
|
|
26 for (j = 0; j < 10; ++j) // mother
|
|
|
27 for (k = 0; k < 10; ++k) // child
|
|
|
28 if ((map[i]&map[k]) && (map[j]&map[k]) && ((map[i]|map[j])&map[k]) == map[k])
|
|
|
29 ret[n++] = j<<16 | i<<8 | k;
|
|
|
30 }
|
|
|
31 }
|
|
|
32 ret[0] = n - 1;
|
|
|
33 return ret;
|
|
|
34 }
|
|
|
35
|
|
|
36
|
|
|
37 int bcf_trio_call(const uint32_t *prep, const bcf1_t *b, int *llr, int64_t *gt)
|
|
|
38 {
|
|
|
39 int i, j, k;
|
|
|
40 const bcf_ginfo_t *PL;
|
|
|
41 uint8_t *gl10;
|
|
|
42 int map[10];
|
|
|
43 if (b->n_smpl != 3) return -1; // not a trio
|
|
|
44 for (i = 0; i < b->n_gi; ++i)
|
|
|
45 if (b->gi[i].fmt == bcf_str2int("PL", 2)) break;
|
|
|
46 if (i == b->n_gi) return -1; // no PL
|
|
|
47 gl10 = alloca(10 * b->n_smpl);
|
|
|
48 if (bcf_gl10(b, gl10) < 0) {
|
|
|
49 if (bcf_gl10_indel(b, gl10) < 0) return -1;
|
|
|
50 }
|
|
|
51 PL = b->gi + i;
|
|
|
52 for (i = 0, k = 0; i < 4; ++i)
|
|
|
53 for (j = i; j < 4; ++j)
|
|
|
54 map[k++] = seq_nt16rev[1<<i|1<<j];
|
|
|
55 for (j = 0; j < 3; ++j) // check if ref hom is the most probable in all members
|
|
|
56 if (((uint8_t*)PL->data)[j * PL->len] != 0) break;
|
|
|
57 if (j < 3) { // we need to go through the complex procedure
|
|
|
58 uint8_t *g[3];
|
|
|
59 int minc = 1<<30, minc_j = -1, minf = 0, gtf = 0, gtc = 0;
|
|
|
60 g[0] = gl10;
|
|
|
61 g[1] = gl10 + 10;
|
|
|
62 g[2] = gl10 + 20;
|
|
|
63 for (j = 1; j <= (int)prep[0]; ++j) { // compute LK with constraint
|
|
|
64 int sum = g[0][prep[j]&0xff] + g[1][prep[j]>>8&0xff] + g[2][prep[j]>>16&0xff];
|
|
|
65 if (sum < minc) minc = sum, minc_j = j;
|
|
|
66 }
|
|
|
67 gtc |= map[prep[minc_j]&0xff]; gtc |= map[prep[minc_j]>>8&0xff]<<8; gtc |= map[prep[minc_j]>>16]<<16;
|
|
|
68 for (j = 0; j < 3; ++j) { // compute LK without constraint
|
|
|
69 int min = 1<<30, min_k = -1;
|
|
|
70 for (k = 0; k < 10; ++k)
|
|
|
71 if (g[j][k] < min) min = g[j][k], min_k = k;
|
|
|
72 gtf |= map[min_k]<<(j*8);
|
|
|
73 minf += min;
|
|
|
74 }
|
|
|
75 *llr = minc - minf; *gt = (int64_t)gtc<<32 | gtf;
|
|
|
76 } else *llr = 0, *gt = -1;
|
|
|
77 return 0;
|
|
|
78 }
|
|
|
79
|
|
|
80 int bcf_pair_call(const bcf1_t *b)
|
|
|
81 {
|
|
|
82 int i, j, k;
|
|
|
83 const bcf_ginfo_t *PL;
|
|
|
84 if (b->n_smpl != 2) return -1; // not a pair
|
|
|
85 for (i = 0; i < b->n_gi; ++i)
|
|
|
86 if (b->gi[i].fmt == bcf_str2int("PL", 2)) break;
|
|
|
87 if (i == b->n_gi) return -1; // no PL
|
|
|
88 PL = b->gi + i;
|
|
|
89 for (j = 0; j < 2; ++j) // check if ref hom is the most probable in all members
|
|
|
90 if (((uint8_t*)PL->data)[j * PL->len] != 0) break;
|
|
|
91 if (j < 2) { // we need to go through the complex procedure
|
|
|
92 uint8_t *g[2];
|
|
|
93 int minc = 1<<30, minf = 0;
|
|
|
94 g[0] = PL->data;
|
|
|
95 g[1] = (uint8_t*)PL->data + PL->len;
|
|
|
96 for (j = 0; j < PL->len; ++j) // compute LK with constraint
|
|
|
97 minc = minc < g[0][j] + g[1][j]? minc : g[0][j] + g[1][j];
|
|
|
98 for (j = 0; j < 2; ++j) { // compute LK without constraint
|
|
|
99 int min = 1<<30;
|
|
|
100 for (k = 0; k < PL->len; ++k)
|
|
|
101 min = min < g[j][k]? min : g[j][k];
|
|
|
102 minf += min;
|
|
|
103 }
|
|
|
104 return minc - minf;
|
|
|
105 } else return 0;
|
|
|
106 }
|
|
|
107
|
|
|
108 int bcf_min_diff(const bcf1_t *b)
|
|
|
109 {
|
|
|
110 int i, min = 1<<30;
|
|
|
111 const bcf_ginfo_t *PL;
|
|
|
112 for (i = 0; i < b->n_gi; ++i)
|
|
|
113 if (b->gi[i].fmt == bcf_str2int("PL", 2)) break;
|
|
|
114 if (i == b->n_gi) return -1; // no PL
|
|
|
115 PL = b->gi + i;
|
|
|
116 for (i = 0; i < b->n_smpl; ++i) {
|
|
|
117 int m1, m2, j;
|
|
|
118 const uint8_t *p = (uint8_t*)PL->data;
|
|
|
119 m1 = m2 = 1<<30;
|
|
|
120 for (j = 0; j < PL->len; ++j) {
|
|
|
121 if ((int)p[j] < m1) m2 = m1, m1 = p[j];
|
|
|
122 else if ((int)p[j] < m2) m2 = p[j];
|
|
|
123 }
|
|
|
124 min = min < m2 - m1? min : m2 - m1;
|
|
|
125 }
|
|
|
126 return min;
|
|
|
127 }
|
