* sequence realignment helps to increase the power.
*/
{ // construct per-sample consensus
- int L = right - left + 1;
- uint32_t *cns;
+ int L = right - left + 1, max_i, max2_i;
+ uint32_t *cns, max, max2;
char *ref0, *r;
ref_sample = calloc(n, sizeof(void*));
cns = calloc(L, 4);
}
}
// determine the consensus
- for (i = 0; i < right - left; ++i)
- r[i] = (cns[i] && (double)(cns[i]&0xffff) / ((cns[i]&0xffff)+(cns[i]>>16&0xffff)) < 0.7)? 15 : ref0[i];
+ for (i = 0; i < right - left; ++i) r[i] = ref0[i];
+ max = max2 = 0; max_i = max2_i = -1;
+ for (i = 0; i < right - left; ++i) {
+ if (cns[i]>>16 >= max>>16) max2 = max, max2_i = max_i, max = cns[i], max_i = i;
+ else if (cns[i]>>16 >= max2>>16) max2 = cns[i], max2_i = i;
+ }
+ if ((double)(max&0xffff) / ((max&0xffff) + (max>>16)) >= 0.7) max_i = -1;
+ if ((double)(max2&0xffff) / ((max2&0xffff) + (max2>>16)) >= 0.7) max2_i = -1;
+ if (max_i >= 0) r[max_i] = 15;
+ if (max2_i >= 0) r[max2_i] = 15;
// for (i = 0; i < right - left; ++i) fputc("=ACMGRSVTWYHKDBN"[(int)r[i]], stderr); fputc('\n', stderr);
}
free(ref0); free(cns);
// align each read to ref2
for (i = 0; i < n_plp[s]; ++i, ++K) {
bam_pileup1_t *p = plp[s] + i;
- int qbeg, qend, tbeg, tend, sc;
+ int qbeg, qend, tbeg, tend, sc, kk;
uint8_t *seq = bam1_seq(p->b);
+ uint32_t *cigar = bam1_cigar(p->b);
+ if (p->b->core.flag&4) continue; // unmapped reads
+ // FIXME: the following loop should be better moved outside; nonetheless, realignment should be much slower anyway.
+ for (kk = 0; kk < p->b->core.n_cigar; ++kk)
+ if ((cigar[kk]&BAM_CIGAR_MASK) == BAM_CREF_SKIP) break;
+ if (kk < p->b->core.n_cigar) continue;
// FIXME: the following skips soft clips, but using them may be more sensitive.
// determine the start and end of sequences for alignment
qbeg = tpos2qpos(&p->b->core, bam1_cigar(p->b), left, 0, &tbeg);
indelQ2 = tmp > 111? 0 : (int)((1. - tmp/111.) * indelQ2 + .499);
// pick the smaller between indelQ1 and indelQ2
indelQ = indelQ1 < indelQ2? indelQ1 : indelQ2;
- p->aux = (sc[0]&0x3f)<<16 | seqQ<<8 | indelQ;
+ p->aux = (sc[0]&0x3f)<<16 | seqQ<<8 | indelQ; // use 22 bits in total
sumq[sc[0]&0x3f] += indelQ < seqQ? indelQ : seqQ;
// fprintf(stderr, "pos=%d read=%d:%d name=%s call=%d q=%d\n", pos, s, i, bam1_qname(p->b), types[sc[0]&0x3f], indelQ);
}