annotate pyPRADA_1.2/tools/samtools-0.1.16/bam_maqcns.c @ 3:f17965495ec9 draft default tip

Uploaded
author siyuan
date Tue, 11 Mar 2014 12:14:01 -0400
parents acc2ca1a3ba4
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
0
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
1 #include <math.h>
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
2 #include <assert.h>
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
3 #include "bam.h"
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
4 #include "bam_maqcns.h"
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
5 #include "ksort.h"
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
6 #include "errmod.h"
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
7 #include "kaln.h"
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
8 KSORT_INIT_GENERIC(uint32_t)
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
9
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
10 #define INDEL_WINDOW_SIZE 50
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
11 #define INDEL_EXT_DEP 0.9
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
12
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
13 typedef struct __bmc_aux_t {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
14 int max;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
15 uint32_t *info;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
16 uint16_t *info16;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
17 errmod_t *em;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
18 } bmc_aux_t;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
19
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
20 typedef struct {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
21 float esum[4], fsum[4];
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
22 uint32_t c[4];
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
23 } glf_call_aux_t;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
24
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
25 /*
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
26 P(<b1,b2>) = \theta \sum_{i=1}^{N-1} 1/i
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
27 P(D|<b1,b2>) = \sum_{k=1}^{N-1} p_k 1/2 [(k/N)^n_2(1-k/N)^n_1 + (k/N)^n1(1-k/N)^n_2]
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
28 p_k = 1/k / \sum_{i=1}^{N-1} 1/i
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
29 */
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
30 static void cal_het(bam_maqcns_t *aa)
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
31 {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
32 int k, n1, n2;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
33 double sum_harmo; // harmonic sum
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
34 double poly_rate;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
35
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
36 free(aa->lhet);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
37 aa->lhet = (double*)calloc(256 * 256, sizeof(double));
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
38 sum_harmo = 0.0;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
39 for (k = 1; k <= aa->n_hap - 1; ++k)
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
40 sum_harmo += 1.0 / k;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
41 for (n1 = 0; n1 < 256; ++n1) {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
42 for (n2 = 0; n2 < 256; ++n2) {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
43 long double sum = 0.0;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
44 double lC = aa->errmod == BAM_ERRMOD_SOAP? 0 : lgamma(n1+n2+1) - lgamma(n1+1) - lgamma(n2+1);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
45 for (k = 1; k <= aa->n_hap - 1; ++k) {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
46 double pk = 1.0 / k / sum_harmo;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
47 double log1 = log((double)k/aa->n_hap);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
48 double log2 = log(1.0 - (double)k/aa->n_hap);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
49 sum += pk * 0.5 * (expl(log1*n2) * expl(log2*n1) + expl(log1*n1) * expl(log2*n2));
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
50 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
51 aa->lhet[n1<<8|n2] = lC + logl(sum);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
52 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
53 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
54 poly_rate = aa->het_rate * sum_harmo;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
55 aa->q_r = -4.343 * log(2.0 * poly_rate / (1.0 - poly_rate));
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
56 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
57
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
58 /** initialize the helper structure */
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
59 static void cal_coef(bam_maqcns_t *aa)
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
60 {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
61 int k, n, q;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
62 long double sum_a[257], b[256], q_c[256], tmp[256], fk2[256];
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
63 double *lC;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
64
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
65 if (aa->errmod == BAM_ERRMOD_MAQ2) return; // no need to do the following
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
66 // aa->lhet will be allocated and initialized
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
67 free(aa->fk); free(aa->coef);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
68 aa->coef = 0;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
69 aa->fk = (double*)calloc(256, sizeof(double));
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
70 aa->fk[0] = fk2[0] = 1.0;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
71 for (n = 1; n != 256; ++n) {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
72 aa->fk[n] = pow(aa->theta, n) * (1.0 - aa->eta) + aa->eta;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
73 fk2[n] = aa->fk[n>>1]; // this is an approximation, assuming reads equally likely come from both strands
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
74 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
75 if (aa->errmod == BAM_ERRMOD_SOAP) return;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
76 aa->coef = (double*)calloc(256*256*64, sizeof(double));
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
77 lC = (double*)calloc(256 * 256, sizeof(double));
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
78 for (n = 1; n != 256; ++n)
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
79 for (k = 1; k <= n; ++k)
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
80 lC[n<<8|k] = lgamma(n+1) - lgamma(k+1) - lgamma(n-k+1);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
81 for (q = 1; q != 64; ++q) {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
82 double e = pow(10.0, -q/10.0);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
83 double le = log(e);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
84 double le1 = log(1.0-e);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
85 for (n = 1; n != 256; ++n) {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
86 double *coef = aa->coef + (q<<16|n<<8);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
87 sum_a[n+1] = 0.0;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
88 for (k = n; k >= 0; --k) { // a_k = \sum_{i=k}^n C^n_k \epsilon^k (1-\epsilon)^{n-k}
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
89 sum_a[k] = sum_a[k+1] + expl(lC[n<<8|k] + k*le + (n-k)*le1);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
90 b[k] = sum_a[k+1] / sum_a[k];
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
91 if (b[k] > 0.99) b[k] = 0.99;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
92 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
93 for (k = 0; k != n; ++k) // log(\bar\beta_{nk}(\bar\epsilon)^{f_k})
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
94 q_c[k] = -4.343 * fk2[k] * logl(b[k] / e);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
95 for (k = 1; k != n; ++k) q_c[k] += q_c[k-1]; // \prod_{i=0}^k c_i
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
96 for (k = 0; k <= n; ++k) { // powl() in 64-bit mode seems broken on my Mac OS X 10.4.9
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
97 tmp[k] = -4.343 * logl(1.0 - expl(fk2[k] * logl(b[k])));
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
98 coef[k] = (k? q_c[k-1] : 0) + tmp[k]; // this is the final c_{nk}
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
99 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
100 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
101 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
102 free(lC);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
103 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
104
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
105 bam_maqcns_t *bam_maqcns_init()
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
106 {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
107 bam_maqcns_t *bm;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
108 bm = (bam_maqcns_t*)calloc(1, sizeof(bam_maqcns_t));
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
109 bm->aux = (bmc_aux_t*)calloc(1, sizeof(bmc_aux_t));
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
110 bm->het_rate = 0.001;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
111 bm->theta = 0.83f;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
112 bm->n_hap = 2;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
113 bm->eta = 0.03;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
114 bm->cap_mapQ = 60;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
115 bm->min_baseQ = 13;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
116 return bm;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
117 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
118
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
119 void bam_maqcns_prepare(bam_maqcns_t *bm)
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
120 {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
121 if (bm->errmod == BAM_ERRMOD_MAQ2) bm->aux->em = errmod_init(1. - bm->theta);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
122 cal_coef(bm); cal_het(bm);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
123 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
124
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
125 void bam_maqcns_destroy(bam_maqcns_t *bm)
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
126 {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
127 if (bm == 0) return;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
128 free(bm->lhet); free(bm->fk); free(bm->coef); free(bm->aux->info); free(bm->aux->info16);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
129 if (bm->aux->em) errmod_destroy(bm->aux->em);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
130 free(bm->aux); free(bm);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
131 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
132
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
133 glf1_t *bam_maqcns_glfgen(int _n, const bam_pileup1_t *pl, uint8_t ref_base, bam_maqcns_t *bm)
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
134 {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
135 glf_call_aux_t *b = 0;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
136 int i, j, k, w[8], c, n;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
137 glf1_t *g = (glf1_t*)calloc(1, sizeof(glf1_t));
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
138 float p[16], min_p = 1e30;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
139 uint64_t rms;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
140
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
141 g->ref_base = ref_base;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
142 if (_n == 0) return g;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
143
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
144 // construct aux array
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
145 if (bm->aux->max < _n) {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
146 bm->aux->max = _n;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
147 kroundup32(bm->aux->max);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
148 bm->aux->info = (uint32_t*)realloc(bm->aux->info, 4 * bm->aux->max);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
149 bm->aux->info16 = (uint16_t*)realloc(bm->aux->info16, 2 * bm->aux->max);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
150 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
151 for (i = n = 0, rms = 0; i < _n; ++i) {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
152 const bam_pileup1_t *p = pl + i;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
153 uint32_t q, x = 0, qq;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
154 uint16_t y = 0;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
155 if (p->is_del || p->is_refskip || (p->b->core.flag&BAM_FUNMAP)) continue;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
156 q = (uint32_t)bam1_qual(p->b)[p->qpos];
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
157 if (q < bm->min_baseQ) continue;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
158 x |= (uint32_t)bam1_strand(p->b) << 18 | q << 8 | p->b->core.qual;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
159 y |= bam1_strand(p->b)<<4;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
160 if (p->b->core.qual < q) q = p->b->core.qual;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
161 c = p->b->core.qual < bm->cap_mapQ? p->b->core.qual : bm->cap_mapQ;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
162 rms += c * c;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
163 x |= q << 24;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
164 y |= q << 5;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
165 qq = bam1_seqi(bam1_seq(p->b), p->qpos);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
166 q = bam_nt16_nt4_table[qq? qq : ref_base];
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
167 if (!p->is_del && !p->is_refskip && q < 4) x |= 1 << 21 | q << 16, y |= q;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
168 bm->aux->info16[n] = y;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
169 bm->aux->info[n++] = x;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
170 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
171 rms = (uint8_t)(sqrt((double)rms / n) + .499);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
172 if (bm->errmod == BAM_ERRMOD_MAQ2) {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
173 errmod_cal(bm->aux->em, n, 4, bm->aux->info16, p);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
174 goto goto_glf;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
175 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
176 ks_introsort(uint32_t, n, bm->aux->info);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
177 // generate esum and fsum
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
178 b = (glf_call_aux_t*)calloc(1, sizeof(glf_call_aux_t));
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
179 for (k = 0; k != 8; ++k) w[k] = 0;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
180 for (j = n - 1; j >= 0; --j) { // calculate esum and fsum
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
181 uint32_t info = bm->aux->info[j];
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
182 if (info>>24 < 4 && (info>>8&0x3f) != 0) info = 4<<24 | (info&0xffffff);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
183 k = info>>16&7;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
184 if (info>>24 > 0) {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
185 b->esum[k&3] += bm->fk[w[k]] * (info>>24);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
186 b->fsum[k&3] += bm->fk[w[k]];
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
187 if (w[k] < 0xff) ++w[k];
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
188 ++b->c[k&3];
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
189 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
190 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
191 // rescale ->c[]
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
192 for (j = c = 0; j != 4; ++j) c += b->c[j];
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
193 if (c > 255) {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
194 for (j = 0; j != 4; ++j) b->c[j] = (int)(254.0 * b->c[j] / c + 0.5);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
195 for (j = c = 0; j != 4; ++j) c += b->c[j];
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
196 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
197 if (bm->errmod == BAM_ERRMOD_MAQ) {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
198 // generate likelihood
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
199 for (j = 0; j != 4; ++j) {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
200 // homozygous
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
201 float tmp1, tmp3;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
202 int tmp2, bar_e;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
203 for (k = 0, tmp1 = tmp3 = 0.0, tmp2 = 0; k != 4; ++k) {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
204 if (j == k) continue;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
205 tmp1 += b->esum[k]; tmp2 += b->c[k]; tmp3 += b->fsum[k];
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
206 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
207 if (tmp2) {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
208 bar_e = (int)(tmp1 / tmp3 + 0.5);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
209 if (bar_e < 4) bar_e = 4; // should not happen
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
210 if (bar_e > 63) bar_e = 63;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
211 p[j<<2|j] = tmp1 + bm->coef[bar_e<<16|c<<8|tmp2];
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
212 } else p[j<<2|j] = 0.0; // all the bases are j
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
213 // heterozygous
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
214 for (k = j + 1; k < 4; ++k) {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
215 for (i = 0, tmp2 = 0, tmp1 = tmp3 = 0.0; i != 4; ++i) {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
216 if (i == j || i == k) continue;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
217 tmp1 += b->esum[i]; tmp2 += b->c[i]; tmp3 += b->fsum[i];
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
218 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
219 if (tmp2) {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
220 bar_e = (int)(tmp1 / tmp3 + 0.5);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
221 if (bar_e < 4) bar_e = 4;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
222 if (bar_e > 63) bar_e = 63;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
223 p[j<<2|k] = p[k<<2|j] = -4.343 * bm->lhet[b->c[j]<<8|b->c[k]] + tmp1 + bm->coef[bar_e<<16|c<<8|tmp2];
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
224 } else p[j<<2|k] = p[k<<2|j] = -4.343 * bm->lhet[b->c[j]<<8|b->c[k]]; // all the bases are either j or k
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
225 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
226 //
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
227 for (k = 0; k != 4; ++k)
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
228 if (p[j<<2|k] < 0.0) p[j<<2|k] = 0.0;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
229 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
230
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
231 { // fix p[k<<2|k]
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
232 float max1, max2, min1, min2;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
233 int max_k, min_k;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
234 max_k = min_k = -1;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
235 max1 = max2 = -1.0; min1 = min2 = 1e30;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
236 for (k = 0; k < 4; ++k) {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
237 if (b->esum[k] > max1) {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
238 max2 = max1; max1 = b->esum[k]; max_k = k;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
239 } else if (b->esum[k] > max2) max2 = b->esum[k];
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
240 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
241 for (k = 0; k < 4; ++k) {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
242 if (p[k<<2|k] < min1) {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
243 min2 = min1; min1 = p[k<<2|k]; min_k = k;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
244 } else if (p[k<<2|k] < min2) min2 = p[k<<2|k];
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
245 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
246 if (max1 > max2 && (min_k != max_k || min1 + 1.0 > min2))
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
247 p[max_k<<2|max_k] = min1 > 1.0? min1 - 1.0 : 0.0;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
248 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
249 } else if (bm->errmod == BAM_ERRMOD_SOAP) { // apply the SOAP model
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
250 // generate likelihood
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
251 for (j = 0; j != 4; ++j) {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
252 float tmp;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
253 // homozygous
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
254 for (k = 0, tmp = 0.0; k != 4; ++k)
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
255 if (j != k) tmp += b->esum[k];
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
256 p[j<<2|j] = tmp;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
257 // heterozygous
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
258 for (k = j + 1; k < 4; ++k) {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
259 for (i = 0, tmp = 0.0; i != 4; ++i)
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
260 if (i != j && i != k) tmp += b->esum[i];
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
261 p[j<<2|k] = p[k<<2|j] = -4.343 * bm->lhet[b->c[j]<<8|b->c[k]] + tmp;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
262 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
263 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
264 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
265
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
266 goto_glf:
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
267 // convert necessary information to glf1_t
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
268 g->ref_base = ref_base; g->max_mapQ = rms;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
269 g->depth = n > 16777215? 16777215 : n;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
270 for (j = 0; j != 4; ++j)
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
271 for (k = j; k < 4; ++k)
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
272 if (p[j<<2|k] < min_p) min_p = p[j<<2|k];
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
273 g->min_lk = min_p > 255.0? 255 : (int)(min_p + 0.5);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
274 for (j = c = 0; j != 4; ++j)
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
275 for (k = j; k < 4; ++k)
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
276 g->lk[c++] = p[j<<2|k]-min_p > 255.0? 255 : (int)(p[j<<2|k]-min_p + 0.5);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
277
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
278 free(b);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
279 return g;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
280 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
281
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
282 uint32_t glf2cns(const glf1_t *g, int q_r)
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
283 {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
284 int i, j, k, p[10], ref4;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
285 uint32_t x = 0;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
286 ref4 = bam_nt16_nt4_table[g->ref_base];
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
287 for (i = k = 0; i < 4; ++i)
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
288 for (j = i; j < 4; ++j) {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
289 int prior = (i == ref4 && j == ref4? 0 : i == ref4 || j == ref4? q_r : q_r + 3);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
290 p[k] = (g->lk[k] + prior)<<4 | i<<2 | j;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
291 ++k;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
292 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
293 for (i = 1; i < 10; ++i) // insertion sort
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
294 for (j = i; j > 0 && p[j] < p[j-1]; --j)
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
295 k = p[j], p[j] = p[j-1], p[j-1] = k;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
296 x = (1u<<(p[0]&3) | 1u<<(p[0]>>2&3)) << 28; // the best genotype
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
297 x |= (uint32_t)g->max_mapQ << 16; // rms mapQ
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
298 x |= ((p[1]>>4) - (p[0]>>4) < 256? (p[1]>>4) - (p[0]>>4) : 255) << 8; // consensus Q
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
299 for (k = 0; k < 10; ++k)
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
300 if ((p[k]&0xf) == (ref4<<2|ref4)) break;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
301 if (k == 10) k = 9;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
302 x |= (p[k]>>4) - (p[0]>>4) < 256? (p[k]>>4) - (p[0]>>4) : 255; // snp Q
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
303 return x;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
304 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
305
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
306 uint32_t bam_maqcns_call(int n, const bam_pileup1_t *pl, bam_maqcns_t *bm)
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
307 {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
308 glf1_t *g;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
309 uint32_t x;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
310 if (n) {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
311 g = bam_maqcns_glfgen(n, pl, 0xf, bm);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
312 x = g->depth == 0? (0xfU<<28 | 0xfU<<24) : glf2cns(g, (int)(bm->q_r + 0.5));
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
313 free(g);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
314 } else x = 0xfU<<28 | 0xfU<<24;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
315 return x;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
316 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
317
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
318 /************** *****************/
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
319
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
320 bam_maqindel_opt_t *bam_maqindel_opt_init()
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
321 {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
322 bam_maqindel_opt_t *mi = (bam_maqindel_opt_t*)calloc(1, sizeof(bam_maqindel_opt_t));
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
323 mi->q_indel = 40;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
324 mi->r_indel = 0.00015;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
325 mi->r_snp = 0.001;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
326 //
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
327 mi->mm_penalty = 3;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
328 mi->indel_err = 4;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
329 mi->ambi_thres = 10;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
330 return mi;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
331 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
332
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
333 void bam_maqindel_ret_destroy(bam_maqindel_ret_t *mir)
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
334 {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
335 if (mir == 0) return;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
336 free(mir->s[0]); free(mir->s[1]); free(mir);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
337 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
338
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
339 int bam_tpos2qpos(const bam1_core_t *c, const uint32_t *cigar, int32_t tpos, int is_left, int32_t *_tpos)
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
340 {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
341 int k, x = c->pos, y = 0, last_y = 0;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
342 *_tpos = c->pos;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
343 for (k = 0; k < c->n_cigar; ++k) {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
344 int op = cigar[k] & BAM_CIGAR_MASK;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
345 int l = cigar[k] >> BAM_CIGAR_SHIFT;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
346 if (op == BAM_CMATCH) {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
347 if (c->pos > tpos) return y;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
348 if (x + l > tpos) {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
349 *_tpos = tpos;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
350 return y + (tpos - x);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
351 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
352 x += l; y += l;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
353 last_y = y;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
354 } else if (op == BAM_CINS || op == BAM_CSOFT_CLIP) y += l;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
355 else if (op == BAM_CDEL || op == BAM_CREF_SKIP) {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
356 if (x + l > tpos) {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
357 *_tpos = is_left? x : x + l;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
358 return y;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
359 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
360 x += l;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
361 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
362 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
363 *_tpos = x;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
364 return last_y;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
365 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
366
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
367 #define MINUS_CONST 0x10000000
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
368
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
369 bam_maqindel_ret_t *bam_maqindel(int n, int pos, const bam_maqindel_opt_t *mi, const bam_pileup1_t *pl, const char *ref,
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
370 int _n_types, int *_types)
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
371 {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
372 int i, j, n_types, *types, left, right, max_rd_len = 0;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
373 bam_maqindel_ret_t *ret = 0;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
374 // if there is no proposed indel, check if there is an indel from the alignment
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
375 if (_n_types == 0) {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
376 for (i = 0; i < n; ++i) {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
377 const bam_pileup1_t *p = pl + i;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
378 if (!(p->b->core.flag&BAM_FUNMAP) && p->indel != 0) break;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
379 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
380 if (i == n) return 0; // no indel
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
381 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
382 { // calculate how many types of indels are available (set n_types and types)
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
383 int m;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
384 uint32_t *aux;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
385 aux = (uint32_t*)calloc(n + _n_types + 1, 4);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
386 m = 0;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
387 aux[m++] = MINUS_CONST; // zero indel is always a type
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
388 for (i = 0; i < n; ++i) {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
389 const bam_pileup1_t *p = pl + i;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
390 if (!(p->b->core.flag&BAM_FUNMAP) && p->indel != 0)
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
391 aux[m++] = MINUS_CONST + p->indel;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
392 j = bam_cigar2qlen(&p->b->core, bam1_cigar(p->b));
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
393 if (j > max_rd_len) max_rd_len = j;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
394 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
395 if (_n_types) // then also add this to aux[]
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
396 for (i = 0; i < _n_types; ++i)
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
397 if (_types[i]) aux[m++] = MINUS_CONST + _types[i];
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
398 ks_introsort(uint32_t, m, aux);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
399 // squeeze out identical types
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
400 for (i = 1, n_types = 1; i < m; ++i)
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
401 if (aux[i] != aux[i-1]) ++n_types;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
402 types = (int*)calloc(n_types, sizeof(int));
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
403 j = 0;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
404 types[j++] = aux[0] - MINUS_CONST;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
405 for (i = 1; i < m; ++i) {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
406 if (aux[i] != aux[i-1])
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
407 types[j++] = aux[i] - MINUS_CONST;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
408 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
409 free(aux);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
410 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
411 { // calculate left and right boundary
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
412 left = pos > INDEL_WINDOW_SIZE? pos - INDEL_WINDOW_SIZE : 0;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
413 right = pos + INDEL_WINDOW_SIZE;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
414 if (types[0] < 0) right -= types[0];
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
415 // in case the alignments stand out the reference
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
416 for (i = pos; i < right; ++i)
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
417 if (ref[i] == 0) break;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
418 right = i;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
419 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
420 { // the core part
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
421 char *ref2, *rs, *inscns = 0;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
422 int qr_snp, k, l, *score, *pscore, max_ins = types[n_types-1];
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
423 qr_snp = (int)(-4.343 * log(mi->r_snp) + .499);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
424 if (max_ins > 0) { // get the consensus of inserted sequences
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
425 int *inscns_aux = (int*)calloc(4 * n_types * max_ins, sizeof(int));
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
426 // count occurrences
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
427 for (i = 0; i < n_types; ++i) {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
428 if (types[i] <= 0) continue; // not insertion
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
429 for (j = 0; j < n; ++j) {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
430 const bam_pileup1_t *p = pl + j;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
431 if (!(p->b->core.flag&BAM_FUNMAP) && p->indel == types[i]) {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
432 for (k = 1; k <= p->indel; ++k) {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
433 int c = bam_nt16_nt4_table[bam1_seqi(bam1_seq(p->b), p->qpos + k)];
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
434 if (c < 4) ++inscns_aux[i*max_ins*4 + (k-1)*4 + c];
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
435 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
436 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
437 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
438 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
439 // construct the consensus of inserted sequence
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
440 inscns = (char*)calloc(n_types * max_ins, sizeof(char));
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
441 for (i = 0; i < n_types; ++i) {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
442 for (j = 0; j < types[i]; ++j) {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
443 int max = 0, max_k = -1, *ia = inscns_aux + i*max_ins*4 + j*4;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
444 for (k = 0; k < 4; ++k) {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
445 if (ia[k] > max) {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
446 max = ia[k];
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
447 max_k = k;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
448 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
449 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
450 inscns[i*max_ins + j] = max? 1<<max_k : 15;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
451 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
452 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
453 free(inscns_aux);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
454 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
455 // calculate score
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
456 ref2 = (char*)calloc(right - left + types[n_types-1] + 2, 1);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
457 rs = (char*)calloc(right - left + max_rd_len + types[n_types-1] + 2, 1);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
458 score = (int*)calloc(n_types * n, sizeof(int));
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
459 pscore = (int*)calloc(n_types * n, sizeof(int));
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
460 for (i = 0; i < n_types; ++i) {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
461 ka_param_t ap = ka_param_blast;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
462 ap.band_width = 2 * types[n_types - 1] + 2;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
463 ap.gap_end_ext = 0;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
464 // write ref2
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
465 for (k = 0, j = left; j <= pos; ++j)
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
466 ref2[k++] = bam_nt16_nt4_table[bam_nt16_table[(int)ref[j]]];
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
467 if (types[i] <= 0) j += -types[i];
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
468 else for (l = 0; l < types[i]; ++l)
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
469 ref2[k++] = bam_nt16_nt4_table[(int)inscns[i*max_ins + l]];
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
470 if (types[0] < 0) { // mask deleted sequences
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
471 int jj, tmp = types[i] >= 0? -types[0] : -types[0] + types[i];
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
472 for (jj = 0; jj < tmp && j < right && ref[j]; ++jj, ++j)
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
473 ref2[k++] = 4;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
474 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
475 for (; j < right && ref[j]; ++j)
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
476 ref2[k++] = bam_nt16_nt4_table[bam_nt16_table[(int)ref[j]]];
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
477 if (j < right) right = j;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
478 // calculate score for each read
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
479 for (j = 0; j < n; ++j) {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
480 const bam_pileup1_t *p = pl + j;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
481 int qbeg, qend, tbeg, tend;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
482 if (p->b->core.flag & BAM_FUNMAP) continue;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
483 qbeg = bam_tpos2qpos(&p->b->core, bam1_cigar(p->b), left, 0, &tbeg);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
484 qend = bam_tpos2qpos(&p->b->core, bam1_cigar(p->b), right, 1, &tend);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
485 assert(tbeg >= left);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
486 for (l = qbeg; l < qend; ++l)
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
487 rs[l - qbeg] = bam_nt16_nt4_table[bam1_seqi(bam1_seq(p->b), l)];
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
488 {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
489 int x, y, n_acigar, ps;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
490 uint32_t *acigar;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
491 ps = 0;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
492 if (tend - tbeg + types[i] <= 0) {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
493 score[i*n+j] = -(1<<20);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
494 pscore[i*n+j] = 1<<20;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
495 continue;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
496 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
497 acigar = ka_global_core((uint8_t*)ref2 + tbeg - left, tend - tbeg + types[i], (uint8_t*)rs, qend - qbeg, &ap, &score[i*n+j], &n_acigar);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
498 x = tbeg - left; y = 0;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
499 for (l = 0; l < n_acigar; ++l) {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
500 int op = acigar[l]&0xf;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
501 int len = acigar[l]>>4;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
502 if (op == BAM_CMATCH) {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
503 int k;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
504 for (k = 0; k < len; ++k)
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
505 if (ref2[x+k] != rs[y+k] && ref2[x+k] < 4)
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
506 ps += bam1_qual(p->b)[y+k] < qr_snp? bam1_qual(p->b)[y+k] : qr_snp;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
507 x += len; y += len;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
508 } else if (op == BAM_CINS || op == BAM_CSOFT_CLIP) {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
509 if (op == BAM_CINS && l > 0 && l < n_acigar - 1) ps += mi->q_indel * len;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
510 y += len;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
511 } else if (op == BAM_CDEL) {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
512 if (l > 0 && l < n_acigar - 1) ps += mi->q_indel * len;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
513 x += len;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
514 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
515 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
516 pscore[i*n+j] = ps;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
517 /*if (1) { // for debugging only
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
518 fprintf(stderr, "id=%d, pos=%d, type=%d, j=%d, score=%d, psore=%d, %d, %d, %d, %d, %d, ",
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
519 j, pos+1, types[i], j, score[i*n+j], pscore[i*n+j], tbeg, tend, qbeg, qend, mi->q_indel);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
520 for (l = 0; l < n_acigar; ++l) fprintf(stderr, "%d%c", acigar[l]>>4, "MIDS"[acigar[l]&0xf]);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
521 fprintf(stderr, "\n");
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
522 for (l = 0; l < tend - tbeg + types[i]; ++l) fputc("ACGTN"[ref2[l+tbeg-left]], stderr);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
523 fputc('\n', stderr);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
524 for (l = 0; l < qend - qbeg; ++l) fputc("ACGTN"[rs[l]], stderr);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
525 fputc('\n', stderr);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
526 }*/
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
527 free(acigar);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
528 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
529 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
530 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
531 { // get final result
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
532 int *sum, max1, max2, max1_i, max2_i;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
533 // pick up the best two score
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
534 sum = (int*)calloc(n_types, sizeof(int));
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
535 for (i = 0; i < n_types; ++i)
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
536 for (j = 0; j < n; ++j)
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
537 sum[i] += -pscore[i*n+j];
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
538 max1 = max2 = -0x7fffffff; max1_i = max2_i = -1;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
539 for (i = 0; i < n_types; ++i) {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
540 if (sum[i] > max1) {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
541 max2 = max1; max2_i = max1_i; max1 = sum[i]; max1_i = i;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
542 } else if (sum[i] > max2) {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
543 max2 = sum[i]; max2_i = i;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
544 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
545 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
546 free(sum);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
547 // write ret
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
548 ret = (bam_maqindel_ret_t*)calloc(1, sizeof(bam_maqindel_ret_t));
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
549 ret->indel1 = types[max1_i]; ret->indel2 = types[max2_i];
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
550 ret->s[0] = (char*)calloc(abs(ret->indel1) + 2, 1);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
551 ret->s[1] = (char*)calloc(abs(ret->indel2) + 2, 1);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
552 // write indel sequence
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
553 if (ret->indel1 > 0) {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
554 ret->s[0][0] = '+';
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
555 for (k = 0; k < ret->indel1; ++k)
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
556 ret->s[0][k+1] = bam_nt16_rev_table[(int)inscns[max1_i*max_ins + k]];
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
557 } else if (ret->indel1 < 0) {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
558 ret->s[0][0] = '-';
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
559 for (k = 0; k < -ret->indel1 && ref[pos + k + 1]; ++k)
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
560 ret->s[0][k+1] = ref[pos + k + 1];
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
561 } else ret->s[0][0] = '*';
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
562 if (ret->indel2 > 0) {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
563 ret->s[1][0] = '+';
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
564 for (k = 0; k < ret->indel2; ++k)
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
565 ret->s[1][k+1] = bam_nt16_rev_table[(int)inscns[max2_i*max_ins + k]];
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
566 } else if (ret->indel2 < 0) {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
567 ret->s[1][0] = '-';
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
568 for (k = 0; k < -ret->indel2 && ref[pos + k + 1]; ++k)
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
569 ret->s[1][k+1] = ref[pos + k + 1];
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
570 } else ret->s[1][0] = '*';
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
571 // write count
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
572 for (i = 0; i < n; ++i) {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
573 const bam_pileup1_t *p = pl + i;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
574 if (p->indel == ret->indel1) ++ret->cnt1;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
575 else if (p->indel == ret->indel2) ++ret->cnt2;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
576 else ++ret->cnt_anti;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
577 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
578 { // write gl[]
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
579 int tmp, seq_err = 0;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
580 double x = 1.0;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
581 tmp = max1_i - max2_i;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
582 if (tmp < 0) tmp = -tmp;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
583 for (j = 0; j < tmp + 1; ++j) x *= INDEL_EXT_DEP;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
584 seq_err = mi->q_indel * (1.0 - x) / (1.0 - INDEL_EXT_DEP);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
585 ret->gl[0] = ret->gl[1] = 0;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
586 for (j = 0; j < n; ++j) {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
587 int s1 = pscore[max1_i*n + j], s2 = pscore[max2_i*n + j];
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
588 //fprintf(stderr, "id=%d, %d, %d, %d, %d, %d\n", j, pl[j].b->core.pos+1, types[max1_i], types[max2_i], s1, s2);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
589 if (s1 > s2) ret->gl[0] += s1 - s2 < seq_err? s1 - s2 : seq_err;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
590 else ret->gl[1] += s2 - s1 < seq_err? s2 - s1 : seq_err;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
591 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
592 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
593 // write cnt_ref and cnt_ambi
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
594 if (max1_i != 0 && max2_i != 0) {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
595 for (j = 0; j < n; ++j) {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
596 int diff1 = score[j] - score[max1_i * n + j];
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
597 int diff2 = score[j] - score[max2_i * n + j];
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
598 if (diff1 > 0 && diff2 > 0) ++ret->cnt_ref;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
599 else if (diff1 == 0 || diff2 == 0) ++ret->cnt_ambi;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
600 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
601 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
602 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
603 free(score); free(pscore); free(ref2); free(rs); free(inscns);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
604 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
605 { // call genotype
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
606 int q[3], qr_indel = (int)(-4.343 * log(mi->r_indel) + 0.5);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
607 int min1, min2, min1_i;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
608 q[0] = ret->gl[0] + (ret->s[0][0] != '*'? 0 : 0) * qr_indel;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
609 q[1] = ret->gl[1] + (ret->s[1][0] != '*'? 0 : 0) * qr_indel;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
610 q[2] = n * 3 + (ret->s[0][0] == '*' || ret->s[1][0] == '*'? 1 : 1) * qr_indel;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
611 min1 = min2 = 0x7fffffff; min1_i = -1;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
612 for (i = 0; i < 3; ++i) {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
613 if (q[i] < min1) {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
614 min2 = min1; min1 = q[i]; min1_i = i;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
615 } else if (q[i] < min2) min2 = q[i];
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
616 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
617 ret->gt = min1_i;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
618 ret->q_cns = min2 - min1;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
619 // set q_ref
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
620 if (ret->gt < 2) ret->q_ref = (ret->s[ret->gt][0] == '*')? 0 : q[1-ret->gt] - q[ret->gt] - qr_indel - 3;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
621 else ret->q_ref = (ret->s[0][0] == '*')? q[0] - q[2] : q[1] - q[2];
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
622 if (ret->q_ref < 0) ret->q_ref = 0;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
623 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
624 free(types);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
625 return ret;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
626 }