// collect read positions for ReadPosBias
int len, pos = get_position(p, &len);
- int epos = (double)pos/len * bca->npos;
+ int epos = (double)pos/(len+1) * bca->npos;
if ( bam1_seqi(bam1_seq(p->b),p->qpos) == ref_base )
bca->ref_pos[epos]++;
else
nref += bca->ref_pos[i];
nalt += bca->alt_pos[i];
U += nref*bca->alt_pos[i];
+ bca->ref_pos[i] = 0;
+ bca->alt_pos[i] = 0;
}
+#if 0
+//todo
double var = 0, avg = (double)(nref+nalt)/bca->npos;
for (i=0; i<bca->npos; i++)
{
call->read_pos.avg = avg;
call->read_pos.var = sqrt(var/bca->npos);
call->read_pos.dp = nref+nalt;
+#endif
if ( !nref || !nalt )
{
call->read_pos_bias = -1;
if (i) kputc(',', &s);
kputw(bc->anno[i], &s);
}
- ksprintf(&s,";RPS=%d,%f,%f", bc->read_pos.dp,bc->read_pos.avg,bc->read_pos.var);
+ //ksprintf(&s,";RPS=%d,%f,%f", bc->read_pos.dp,bc->read_pos.avg,bc->read_pos.var);
ksprintf(&s,";QS=%f,%f,%f,%f", bc->qsum[0],bc->qsum[1],bc->qsum[2],bc->qsum[3]);
if (bc->vdb != -1)
ksprintf(&s, ";VDB=%e", bc->vdb);