|
0
|
1 #include <math.h>
|
|
|
2 #include <stdint.h>
|
|
|
3 #include <assert.h>
|
|
|
4 #include "bam.h"
|
|
|
5 #include "kstring.h"
|
|
|
6 #include "bam2bcf.h"
|
|
|
7 #include "errmod.h"
|
|
|
8 #include "bcftools/bcf.h"
|
|
|
9
|
|
|
10 extern void ks_introsort_uint32_t(size_t n, uint32_t a[]);
|
|
|
11
|
|
|
12 #define CALL_ETA 0.03f
|
|
|
13 #define CALL_MAX 256
|
|
|
14 #define CALL_DEFTHETA 0.83f
|
|
|
15 #define DEF_MAPQ 20
|
|
|
16
|
|
|
17 #define CAP_DIST 25
|
|
|
18
|
|
|
19 bcf_callaux_t *bcf_call_init(double theta, int min_baseQ)
|
|
|
20 {
|
|
|
21 bcf_callaux_t *bca;
|
|
|
22 if (theta <= 0.) theta = CALL_DEFTHETA;
|
|
|
23 bca = calloc(1, sizeof(bcf_callaux_t));
|
|
|
24 bca->capQ = 60;
|
|
|
25 bca->openQ = 40; bca->extQ = 20; bca->tandemQ = 100;
|
|
|
26 bca->min_baseQ = min_baseQ;
|
|
|
27 bca->e = errmod_init(1. - theta);
|
|
|
28 bca->min_frac = 0.002;
|
|
|
29 bca->min_support = 1;
|
|
|
30 bca->per_sample_flt = 0;
|
|
|
31 bca->npos = 100;
|
|
|
32 bca->ref_pos = calloc(bca->npos, sizeof(int));
|
|
|
33 bca->alt_pos = calloc(bca->npos, sizeof(int));
|
|
|
34 return bca;
|
|
|
35 }
|
|
|
36
|
|
|
37
|
|
|
38 static int get_position(const bam_pileup1_t *p, int *len)
|
|
|
39 {
|
|
|
40 int icig, n_tot_bases = 0, iread = 0, edist = p->qpos + 1;
|
|
|
41 for (icig=0; icig<p->b->core.n_cigar; icig++)
|
|
|
42 {
|
|
|
43 // Conversion from uint32_t to MIDNSHP
|
|
|
44 // 0123456
|
|
|
45 // MIDNSHP
|
|
|
46 int cig = bam1_cigar(p->b)[icig] & BAM_CIGAR_MASK;
|
|
|
47 int ncig = bam1_cigar(p->b)[icig] >> BAM_CIGAR_SHIFT;
|
|
|
48 if ( cig==0 )
|
|
|
49 {
|
|
|
50 n_tot_bases += ncig;
|
|
|
51 iread += ncig;
|
|
|
52 }
|
|
|
53 else if ( cig==1 )
|
|
|
54 {
|
|
|
55 n_tot_bases += ncig;
|
|
|
56 iread += ncig;
|
|
|
57 }
|
|
|
58 else if ( cig==4 )
|
|
|
59 {
|
|
|
60 iread += ncig;
|
|
|
61 if ( iread<=p->qpos ) edist -= ncig;
|
|
|
62 }
|
|
|
63 }
|
|
|
64 *len = n_tot_bases;
|
|
|
65 return edist;
|
|
|
66 }
|
|
|
67
|
|
|
68 void bcf_call_destroy(bcf_callaux_t *bca)
|
|
|
69 {
|
|
|
70 if (bca == 0) return;
|
|
|
71 errmod_destroy(bca->e);
|
|
|
72 if (bca->npos) { free(bca->ref_pos); free(bca->alt_pos); bca->npos = 0; }
|
|
|
73 free(bca->bases); free(bca->inscns); free(bca);
|
|
|
74 }
|
|
|
75 /* ref_base is the 4-bit representation of the reference base. It is
|
|
|
76 * negative if we are looking at an indel. */
|
|
|
77 int bcf_call_glfgen(int _n, const bam_pileup1_t *pl, int ref_base, bcf_callaux_t *bca, bcf_callret1_t *r)
|
|
|
78 {
|
|
|
79 int i, n, ref4, is_indel, ori_depth = 0;
|
|
|
80 memset(r, 0, sizeof(bcf_callret1_t));
|
|
|
81 if (ref_base >= 0) {
|
|
|
82 ref4 = bam_nt16_nt4_table[ref_base];
|
|
|
83 is_indel = 0;
|
|
|
84 } else ref4 = 4, is_indel = 1;
|
|
|
85 if (_n == 0) return -1;
|
|
|
86 // enlarge the bases array if necessary
|
|
|
87 if (bca->max_bases < _n) {
|
|
|
88 bca->max_bases = _n;
|
|
|
89 kroundup32(bca->max_bases);
|
|
|
90 bca->bases = (uint16_t*)realloc(bca->bases, 2 * bca->max_bases);
|
|
|
91 }
|
|
|
92 // fill the bases array
|
|
|
93 for (i = n = r->n_supp = 0; i < _n; ++i) {
|
|
|
94 const bam_pileup1_t *p = pl + i;
|
|
|
95 int q, b, mapQ, baseQ, is_diff, min_dist, seqQ;
|
|
|
96 // set base
|
|
|
97 if (p->is_del || p->is_refskip || (p->b->core.flag&BAM_FUNMAP)) continue;
|
|
|
98 ++ori_depth;
|
|
|
99 baseQ = q = is_indel? p->aux&0xff : (int)bam1_qual(p->b)[p->qpos]; // base/indel quality
|
|
|
100 seqQ = is_indel? (p->aux>>8&0xff) : 99;
|
|
|
101 if (q < bca->min_baseQ) continue;
|
|
|
102 if (q > seqQ) q = seqQ;
|
|
|
103 mapQ = p->b->core.qual < 255? p->b->core.qual : DEF_MAPQ; // special case for mapQ==255
|
|
|
104 mapQ = mapQ < bca->capQ? mapQ : bca->capQ;
|
|
|
105 if (q > mapQ) q = mapQ;
|
|
|
106 if (q > 63) q = 63;
|
|
|
107 if (q < 4) q = 4;
|
|
|
108 if (!is_indel) {
|
|
|
109 b = bam1_seqi(bam1_seq(p->b), p->qpos); // base
|
|
|
110 b = bam_nt16_nt4_table[b? b : ref_base]; // b is the 2-bit base
|
|
|
111 is_diff = (ref4 < 4 && b == ref4)? 0 : 1;
|
|
|
112 } else {
|
|
|
113 b = p->aux>>16&0x3f;
|
|
|
114 is_diff = (b != 0);
|
|
|
115 }
|
|
|
116 if (is_diff) ++r->n_supp;
|
|
|
117 bca->bases[n++] = q<<5 | (int)bam1_strand(p->b)<<4 | b;
|
|
|
118 // collect annotations
|
|
|
119 if (b < 4) r->qsum[b] += q;
|
|
|
120 ++r->anno[0<<2|is_diff<<1|bam1_strand(p->b)];
|
|
|
121 min_dist = p->b->core.l_qseq - 1 - p->qpos;
|
|
|
122 if (min_dist > p->qpos) min_dist = p->qpos;
|
|
|
123 if (min_dist > CAP_DIST) min_dist = CAP_DIST;
|
|
|
124 r->anno[1<<2|is_diff<<1|0] += baseQ;
|
|
|
125 r->anno[1<<2|is_diff<<1|1] += baseQ * baseQ;
|
|
|
126 r->anno[2<<2|is_diff<<1|0] += mapQ;
|
|
|
127 r->anno[2<<2|is_diff<<1|1] += mapQ * mapQ;
|
|
|
128 r->anno[3<<2|is_diff<<1|0] += min_dist;
|
|
|
129 r->anno[3<<2|is_diff<<1|1] += min_dist * min_dist;
|
|
|
130
|
|
|
131 // collect read positions for ReadPosBias
|
|
|
132 int len, pos = get_position(p, &len);
|
|
|
133 int epos = (double)pos/(len+1) * bca->npos;
|
|
|
134 if ( bam1_seqi(bam1_seq(p->b),p->qpos) == ref_base )
|
|
|
135 bca->ref_pos[epos]++;
|
|
|
136 else
|
|
|
137 bca->alt_pos[epos]++;
|
|
|
138 }
|
|
|
139 r->depth = n; r->ori_depth = ori_depth;
|
|
|
140 // glfgen
|
|
|
141 errmod_cal(bca->e, n, 5, bca->bases, r->p);
|
|
|
142 return r->depth;
|
|
|
143 }
|
|
|
144
|
|
|
145 double mann_whitney_1947(int n, int m, int U)
|
|
|
146 {
|
|
|
147 if (U<0) return 0;
|
|
|
148 if (n==0||m==0) return U==0 ? 1 : 0;
|
|
|
149 return (double)n/(n+m)*mann_whitney_1947(n-1,m,U-m) + (double)m/(n+m)*mann_whitney_1947(n,m-1,U);
|
|
|
150 }
|
|
|
151
|
|
|
152 void calc_ReadPosBias(bcf_callaux_t *bca, bcf_call_t *call)
|
|
|
153 {
|
|
|
154 int i, nref = 0, nalt = 0;
|
|
|
155 unsigned long int U = 0;
|
|
|
156 for (i=0; i<bca->npos; i++)
|
|
|
157 {
|
|
|
158 nref += bca->ref_pos[i];
|
|
|
159 nalt += bca->alt_pos[i];
|
|
|
160 U += nref*bca->alt_pos[i];
|
|
|
161 bca->ref_pos[i] = 0;
|
|
|
162 bca->alt_pos[i] = 0;
|
|
|
163 }
|
|
|
164 #if 0
|
|
|
165 //todo
|
|
|
166 double var = 0, avg = (double)(nref+nalt)/bca->npos;
|
|
|
167 for (i=0; i<bca->npos; i++)
|
|
|
168 {
|
|
|
169 double ediff = bca->ref_pos[i] + bca->alt_pos[i] - avg;
|
|
|
170 var += ediff*ediff;
|
|
|
171 bca->ref_pos[i] = 0;
|
|
|
172 bca->alt_pos[i] = 0;
|
|
|
173 }
|
|
|
174 call->read_pos.avg = avg;
|
|
|
175 call->read_pos.var = sqrt(var/bca->npos);
|
|
|
176 call->read_pos.dp = nref+nalt;
|
|
|
177 #endif
|
|
|
178 if ( !nref || !nalt )
|
|
|
179 {
|
|
|
180 call->read_pos_bias = -1;
|
|
|
181 return;
|
|
|
182 }
|
|
|
183
|
|
|
184 if ( nref>=8 || nalt>=8 )
|
|
|
185 {
|
|
|
186 // normal approximation
|
|
|
187 double mean = ((double)nref*nalt+1.0)/2.0;
|
|
|
188 double var2 = (double)nref*nalt*(nref+nalt+1.0)/12.0;
|
|
|
189 double z = (U-mean)/sqrt(var2);
|
|
|
190 call->read_pos_bias = z;
|
|
|
191 //fprintf(stderr,"nref=%d nalt=%d U=%ld mean=%e var=%e zval=%e\n", nref,nalt,U,mean,sqrt(var2),call->read_pos_bias);
|
|
|
192 }
|
|
|
193 else
|
|
|
194 {
|
|
|
195 double p = mann_whitney_1947(nalt,nref,U);
|
|
|
196 // biased form claimed by GATK to behave better empirically
|
|
|
197 // double var2 = (1.0+1.0/(nref+nalt+1.0))*(double)nref*nalt*(nref+nalt+1.0)/12.0;
|
|
|
198 double var2 = (double)nref*nalt*(nref+nalt+1.0)/12.0;
|
|
|
199 double z;
|
|
|
200 if ( p >= 1./sqrt(var2*2*M_PI) ) z = 0; // equal to mean
|
|
|
201 else
|
|
|
202 {
|
|
|
203 if ( U >= nref*nalt/2. ) z = sqrt(-2*log(sqrt(var2*2*M_PI)*p));
|
|
|
204 else z = -sqrt(-2*log(sqrt(var2*2*M_PI)*p));
|
|
|
205 }
|
|
|
206 call->read_pos_bias = z;
|
|
|
207 //fprintf(stderr,"nref=%d nalt=%d U=%ld p=%e var2=%e zval=%e\n", nref,nalt,U, p,var2,call->read_pos_bias);
|
|
|
208 }
|
|
|
209 }
|
|
|
210
|
|
|
211 float mean_diff_to_prob(float mdiff, int dp, int readlen)
|
|
|
212 {
|
|
|
213 if ( dp==2 )
|
|
|
214 {
|
|
|
215 if ( mdiff==0 )
|
|
|
216 return (2.0*readlen + 4.0*(readlen-1.0))/((float)readlen*readlen);
|
|
|
217 else
|
|
|
218 return 8.0*(readlen - 4.0*mdiff)/((float)readlen*readlen);
|
|
|
219 }
|
|
|
220
|
|
|
221 // This is crude empirical approximation and is not very accurate for
|
|
|
222 // shorter read lengths (<100bp). There certainly is a room for
|
|
|
223 // improvement.
|
|
|
224 const float mv[24][2] = { {0,0}, {0,0}, {0,0},
|
|
|
225 { 9.108, 4.934}, { 9.999, 3.991}, {10.273, 3.485}, {10.579, 3.160},
|
|
|
226 {10.828, 2.889}, {11.014, 2.703}, {11.028, 2.546}, {11.244, 2.391},
|
|
|
227 {11.231, 2.320}, {11.323, 2.138}, {11.403, 2.123}, {11.394, 1.994},
|
|
|
228 {11.451, 1.928}, {11.445, 1.862}, {11.516, 1.815}, {11.560, 1.761},
|
|
|
229 {11.544, 1.728}, {11.605, 1.674}, {11.592, 1.652}, {11.674, 1.613},
|
|
|
230 {11.641, 1.570} };
|
|
|
231
|
|
|
232 float m, v;
|
|
|
233 if ( dp>=24 )
|
|
|
234 {
|
|
|
235 m = readlen/8.;
|
|
|
236 if (dp>100) dp = 100;
|
|
|
237 v = 1.476/(0.182*pow(dp,0.514));
|
|
|
238 v = v*(readlen/100.);
|
|
|
239 }
|
|
|
240 else
|
|
|
241 {
|
|
|
242 m = mv[dp][0];
|
|
|
243 v = mv[dp][1];
|
|
|
244 m = m*readlen/100.;
|
|
|
245 v = v*readlen/100.;
|
|
|
246 v *= 1.2; // allow more variability
|
|
|
247 }
|
|
|
248 return 1.0/(v*sqrt(2*M_PI)) * exp(-0.5*((mdiff-m)/v)*((mdiff-m)/v));
|
|
|
249 }
|
|
|
250
|
|
|
251 void calc_vdb(bcf_callaux_t *bca, bcf_call_t *call)
|
|
|
252 {
|
|
|
253 int i, dp = 0;
|
|
|
254 float mean_pos = 0, mean_diff = 0;
|
|
|
255 for (i=0; i<bca->npos; i++)
|
|
|
256 {
|
|
|
257 if ( !bca->alt_pos[i] ) continue;
|
|
|
258 dp += bca->alt_pos[i];
|
|
|
259 int j = i<bca->npos/2 ? i : bca->npos - i;
|
|
|
260 mean_pos += bca->alt_pos[i]*j;
|
|
|
261 }
|
|
|
262 if ( dp<2 )
|
|
|
263 {
|
|
|
264 call->vdb = -1;
|
|
|
265 return;
|
|
|
266 }
|
|
|
267 mean_pos /= dp;
|
|
|
268 for (i=0; i<bca->npos; i++)
|
|
|
269 {
|
|
|
270 if ( !bca->alt_pos[i] ) continue;
|
|
|
271 int j = i<bca->npos/2 ? i : bca->npos - i;
|
|
|
272 mean_diff += bca->alt_pos[i] * fabs(j - mean_pos);
|
|
|
273 }
|
|
|
274 mean_diff /= dp;
|
|
|
275 call->vdb = mean_diff_to_prob(mean_diff, dp, bca->npos);
|
|
|
276 }
|
|
|
277
|
|
|
278 /**
|
|
|
279 * bcf_call_combine() - sets the PL array and VDB, RPB annotations, finds the top two alleles
|
|
|
280 * @n: number of samples
|
|
|
281 * @calls: each sample's calls
|
|
|
282 * @bca: auxiliary data structure for holding temporary values
|
|
|
283 * @ref_base: the reference base
|
|
|
284 * @call: filled with the annotations
|
|
|
285 */
|
|
|
286 int bcf_call_combine(int n, const bcf_callret1_t *calls, bcf_callaux_t *bca, int ref_base /*4-bit*/, bcf_call_t *call)
|
|
|
287 {
|
|
|
288 int ref4, i, j, qsum[4];
|
|
|
289 int64_t tmp;
|
|
|
290 if (ref_base >= 0) {
|
|
|
291 call->ori_ref = ref4 = bam_nt16_nt4_table[ref_base];
|
|
|
292 if (ref4 > 4) ref4 = 4;
|
|
|
293 } else call->ori_ref = -1, ref4 = 0;
|
|
|
294 // calculate qsum
|
|
|
295 memset(qsum, 0, 4 * sizeof(int));
|
|
|
296 for (i = 0; i < n; ++i)
|
|
|
297 for (j = 0; j < 4; ++j)
|
|
|
298 qsum[j] += calls[i].qsum[j];
|
|
|
299 int qsum_tot=0;
|
|
|
300 for (j=0; j<4; j++) { qsum_tot += qsum[j]; call->qsum[j] = 0; }
|
|
|
301 for (j = 0; j < 4; ++j) qsum[j] = qsum[j] << 2 | j;
|
|
|
302 // find the top 2 alleles
|
|
|
303 for (i = 1; i < 4; ++i) // insertion sort
|
|
|
304 for (j = i; j > 0 && qsum[j] < qsum[j-1]; --j)
|
|
|
305 tmp = qsum[j], qsum[j] = qsum[j-1], qsum[j-1] = tmp;
|
|
|
306 // set the reference allele and alternative allele(s)
|
|
|
307 for (i = 0; i < 5; ++i) call->a[i] = -1;
|
|
|
308 call->unseen = -1;
|
|
|
309 call->a[0] = ref4;
|
|
|
310 for (i = 3, j = 1; i >= 0; --i) {
|
|
|
311 if ((qsum[i]&3) != ref4) {
|
|
|
312 if (qsum[i]>>2 != 0)
|
|
|
313 {
|
|
|
314 if ( j<4 ) call->qsum[j] = (float)(qsum[i]>>2)/qsum_tot; // ref N can make j>=4
|
|
|
315 call->a[j++] = qsum[i]&3;
|
|
|
316 }
|
|
|
317 else break;
|
|
|
318 }
|
|
|
319 else
|
|
|
320 call->qsum[0] = (float)(qsum[i]>>2)/qsum_tot;
|
|
|
321 }
|
|
|
322 if (ref_base >= 0) { // for SNPs, find the "unseen" base
|
|
|
323 if (((ref4 < 4 && j < 4) || (ref4 == 4 && j < 5)) && i >= 0)
|
|
|
324 call->unseen = j, call->a[j++] = qsum[i]&3;
|
|
|
325 call->n_alleles = j;
|
|
|
326 } else {
|
|
|
327 call->n_alleles = j;
|
|
|
328 if (call->n_alleles == 1) return -1; // no reliable supporting read. stop doing anything
|
|
|
329 }
|
|
|
330 // set the PL array
|
|
|
331 if (call->n < n) {
|
|
|
332 call->n = n;
|
|
|
333 call->PL = realloc(call->PL, 15 * n);
|
|
|
334 }
|
|
|
335 {
|
|
|
336 int x, g[15], z;
|
|
|
337 double sum_min = 0.;
|
|
|
338 x = call->n_alleles * (call->n_alleles + 1) / 2;
|
|
|
339 // get the possible genotypes
|
|
|
340 for (i = z = 0; i < call->n_alleles; ++i)
|
|
|
341 for (j = 0; j <= i; ++j)
|
|
|
342 g[z++] = call->a[j] * 5 + call->a[i];
|
|
|
343 for (i = 0; i < n; ++i) {
|
|
|
344 uint8_t *PL = call->PL + x * i;
|
|
|
345 const bcf_callret1_t *r = calls + i;
|
|
|
346 float min = 1e37;
|
|
|
347 for (j = 0; j < x; ++j)
|
|
|
348 if (min > r->p[g[j]]) min = r->p[g[j]];
|
|
|
349 sum_min += min;
|
|
|
350 for (j = 0; j < x; ++j) {
|
|
|
351 int y;
|
|
|
352 y = (int)(r->p[g[j]] - min + .499);
|
|
|
353 if (y > 255) y = 255;
|
|
|
354 PL[j] = y;
|
|
|
355 }
|
|
|
356 }
|
|
|
357 // if (ref_base < 0) fprintf(stderr, "%d,%d,%f,%d\n", call->n_alleles, x, sum_min, call->unseen);
|
|
|
358 call->shift = (int)(sum_min + .499);
|
|
|
359 }
|
|
|
360 // combine annotations
|
|
|
361 memset(call->anno, 0, 16 * sizeof(int));
|
|
|
362 for (i = call->depth = call->ori_depth = 0, tmp = 0; i < n; ++i) {
|
|
|
363 call->depth += calls[i].depth;
|
|
|
364 call->ori_depth += calls[i].ori_depth;
|
|
|
365 for (j = 0; j < 16; ++j) call->anno[j] += calls[i].anno[j];
|
|
|
366 }
|
|
|
367
|
|
|
368 calc_vdb(bca, call);
|
|
|
369 calc_ReadPosBias(bca, call);
|
|
|
370
|
|
|
371 return 0;
|
|
|
372 }
|
|
|
373
|
|
|
374 int bcf_call2bcf(int tid, int pos, bcf_call_t *bc, bcf1_t *b, bcf_callret1_t *bcr, int fmt_flag,
|
|
|
375 const bcf_callaux_t *bca, const char *ref)
|
|
|
376 {
|
|
|
377 extern double kt_fisher_exact(int n11, int n12, int n21, int n22, double *_left, double *_right, double *two);
|
|
|
378 kstring_t s;
|
|
|
379 int i, j;
|
|
|
380 b->n_smpl = bc->n;
|
|
|
381 b->tid = tid; b->pos = pos; b->qual = 0;
|
|
|
382 s.s = b->str; s.m = b->m_str; s.l = 0;
|
|
|
383 kputc('\0', &s);
|
|
|
384 if (bc->ori_ref < 0) { // an indel
|
|
|
385 // write REF
|
|
|
386 kputc(ref[pos], &s);
|
|
|
387 for (j = 0; j < bca->indelreg; ++j) kputc(ref[pos+1+j], &s);
|
|
|
388 kputc('\0', &s);
|
|
|
389 // write ALT
|
|
|
390 kputc(ref[pos], &s);
|
|
|
391 for (i = 1; i < 4; ++i) {
|
|
|
392 if (bc->a[i] < 0) break;
|
|
|
393 if (i > 1) {
|
|
|
394 kputc(',', &s); kputc(ref[pos], &s);
|
|
|
395 }
|
|
|
396 if (bca->indel_types[bc->a[i]] < 0) { // deletion
|
|
|
397 for (j = -bca->indel_types[bc->a[i]]; j < bca->indelreg; ++j)
|
|
|
398 kputc(ref[pos+1+j], &s);
|
|
|
399 } else { // insertion; cannot be a reference unless a bug
|
|
|
400 char *inscns = &bca->inscns[bc->a[i] * bca->maxins];
|
|
|
401 for (j = 0; j < bca->indel_types[bc->a[i]]; ++j)
|
|
|
402 kputc("ACGTN"[(int)inscns[j]], &s);
|
|
|
403 for (j = 0; j < bca->indelreg; ++j) kputc(ref[pos+1+j], &s);
|
|
|
404 }
|
|
|
405 }
|
|
|
406 kputc('\0', &s);
|
|
|
407 } else { // a SNP
|
|
|
408 kputc("ACGTN"[bc->ori_ref], &s); kputc('\0', &s);
|
|
|
409 for (i = 1; i < 5; ++i) {
|
|
|
410 if (bc->a[i] < 0) break;
|
|
|
411 if (i > 1) kputc(',', &s);
|
|
|
412 kputc(bc->unseen == i? 'X' : "ACGT"[bc->a[i]], &s);
|
|
|
413 }
|
|
|
414 kputc('\0', &s);
|
|
|
415 }
|
|
|
416 kputc('\0', &s);
|
|
|
417 // INFO
|
|
|
418 if (bc->ori_ref < 0) ksprintf(&s,"INDEL;IS=%d,%f;", bca->max_support, bca->max_frac);
|
|
|
419 kputs("DP=", &s); kputw(bc->ori_depth, &s); kputs(";I16=", &s);
|
|
|
420 for (i = 0; i < 16; ++i) {
|
|
|
421 if (i) kputc(',', &s);
|
|
|
422 kputw(bc->anno[i], &s);
|
|
|
423 }
|
|
|
424 //ksprintf(&s,";RPS=%d,%f,%f", bc->read_pos.dp,bc->read_pos.avg,bc->read_pos.var);
|
|
|
425 ksprintf(&s,";QS=%f,%f,%f,%f", bc->qsum[0],bc->qsum[1],bc->qsum[2],bc->qsum[3]);
|
|
|
426 if (bc->vdb != -1)
|
|
|
427 ksprintf(&s, ";VDB=%e", bc->vdb);
|
|
|
428 if (bc->read_pos_bias != -1 )
|
|
|
429 ksprintf(&s, ";RPB=%e", bc->read_pos_bias);
|
|
|
430 kputc('\0', &s);
|
|
|
431 // FMT
|
|
|
432 kputs("PL", &s);
|
|
|
433 if (bcr && fmt_flag) {
|
|
|
434 if (fmt_flag & B2B_FMT_DP) kputs(":DP", &s);
|
|
|
435 if (fmt_flag & B2B_FMT_DV) kputs(":DV", &s);
|
|
|
436 if (fmt_flag & B2B_FMT_SP) kputs(":SP", &s);
|
|
|
437 }
|
|
|
438 kputc('\0', &s);
|
|
|
439 b->m_str = s.m; b->str = s.s; b->l_str = s.l;
|
|
|
440 bcf_sync(b);
|
|
|
441 memcpy(b->gi[0].data, bc->PL, b->gi[0].len * bc->n);
|
|
|
442 if (bcr && fmt_flag) {
|
|
|
443 uint16_t *dp = (fmt_flag & B2B_FMT_DP)? b->gi[1].data : 0;
|
|
|
444 uint16_t *dv = (fmt_flag & B2B_FMT_DV)? b->gi[1 + ((fmt_flag & B2B_FMT_DP) != 0)].data : 0;
|
|
|
445 int32_t *sp = (fmt_flag & B2B_FMT_SP)? b->gi[1 + ((fmt_flag & B2B_FMT_DP) != 0) + ((fmt_flag & B2B_FMT_DV) != 0)].data : 0;
|
|
|
446 for (i = 0; i < bc->n; ++i) {
|
|
|
447 bcf_callret1_t *p = bcr + i;
|
|
|
448 if (dp) dp[i] = p->depth < 0xffff? p->depth : 0xffff;
|
|
|
449 if (dv) dv[i] = p->n_supp < 0xffff? p->n_supp : 0xffff;
|
|
|
450 if (sp) {
|
|
|
451 if (p->anno[0] + p->anno[1] < 2 || p->anno[2] + p->anno[3] < 2
|
|
|
452 || p->anno[0] + p->anno[2] < 2 || p->anno[1] + p->anno[3] < 2)
|
|
|
453 {
|
|
|
454 sp[i] = 0;
|
|
|
455 } else {
|
|
|
456 double left, right, two;
|
|
|
457 int x;
|
|
|
458 kt_fisher_exact(p->anno[0], p->anno[1], p->anno[2], p->anno[3], &left, &right, &two);
|
|
|
459 x = (int)(-4.343 * log(two) + .499);
|
|
|
460 if (x > 255) x = 255;
|
|
|
461 sp[i] = x;
|
|
|
462 }
|
|
|
463 }
|
|
|
464 }
|
|
|
465 }
|
|
|
466 return 0;
|
|
|
467 }
|
