if (theta <= 0.) theta = CALL_DEFTHETA;
bca = calloc(1, sizeof(bcf_callaux_t));
bca->capQ = 60;
- bca->openQ = 40; bca->extQ = 20; bca->tandemQ = 100;
+ bca->openQ = 40; bca->extQ = 20; bca->tandemQ = 80;
bca->min_baseQ = min_baseQ;
bca->e = errmod_init(1. - theta);
return bca;
errmod_destroy(bca->e);
free(bca->bases); free(bca);
}
-
-int bcf_call_glfgen(int _n, const bam_pileup1_t *pl, int ref_base /*4-bit*/, bcf_callaux_t *bca, bcf_callret1_t *r)
+/* ref_base is the 4-bit representation of the reference base. It is
+ * negative if we are looking at an indel. */
+int bcf_call_glfgen(int _n, const bam_pileup1_t *pl, int ref_base, bcf_callaux_t *bca, bcf_callret1_t *r)
{
- int i, n, ref4;
+ int i, n, ref4, is_indel;
memset(r, 0, sizeof(bcf_callret1_t));
- ref4 = bam_nt16_nt4_table[ref_base];
+ if (ref_base >= 0) {
+ ref4 = bam_nt16_nt4_table[ref_base];
+ is_indel = 0;
+ } else ref4 = 4, is_indel = 1;
if (_n == 0) return -1;
-
// enlarge the bases array if necessary
if (bca->max_bases < _n) {
bca->max_bases = _n;
const bam_pileup1_t *p = pl + i;
int q, b, mapQ, baseQ, is_diff, min_dist;
// set base
- if (p->is_del || (p->b->core.flag&BAM_FUNMAP)) continue; // skip unmapped reads and deleted bases
- baseQ = q = (int)bam1_qual(p->b)[p->qpos]; // base quality
+ if (p->is_del || p->is_refskip || (p->b->core.flag&BAM_FUNMAP)) continue;
+ baseQ = q = is_indel? p->aux&0xff : (int)bam1_qual(p->b)[p->qpos]; // base/indel quality
if (q < bca->min_baseQ) continue;
mapQ = p->b->core.qual < bca->capQ? p->b->core.qual : bca->capQ;
if (q > mapQ) q = mapQ;
if (q > 63) q = 63;
if (q < 4) q = 4;
- b = bam1_seqi(bam1_seq(p->b), p->qpos); // base
- b = bam_nt16_nt4_table[b? b : ref_base]; // b is the 2-bit base
+ if (!is_indel) {
+ b = bam1_seqi(bam1_seq(p->b), p->qpos); // base
+ b = bam_nt16_nt4_table[b? b : ref_base]; // b is the 2-bit base
+ is_diff = (ref4 < 4 && b == ref4)? 0 : 1;
+ } else {
+ b = p->aux>>8&0x3f;
+ is_diff = (b != 0);
+ }
bca->bases[n++] = q<<5 | (int)bam1_strand(p->b)<<4 | b;
// collect annotations
r->qsum[b] += q;
- is_diff = (ref4 < 4 && b == ref4)? 0 : 1;
++r->anno[0<<2|is_diff<<1|bam1_strand(p->b)];
min_dist = p->b->core.l_qseq - 1 - p->qpos;
if (min_dist > p->qpos) min_dist = p->qpos;
{
int ref4, i, j, qsum[4];
int64_t tmp;
- call->ori_ref = ref4 = bam_nt16_nt4_table[ref_base];
- if (ref4 > 4) ref4 = 4;
+ if (ref_base >= 0) {
+ call->ori_ref = ref4 = bam_nt16_nt4_table[ref_base];
+ if (ref4 > 4) ref4 = 4;
+ } else call->ori_ref = -1, ref4 = 0;
// calculate qsum
memset(qsum, 0, 4 * sizeof(int));
for (i = 0; i < n; ++i)
else break;
}
}
- if (((ref4 < 4 && j < 4) || (ref4 == 4 && j < 5)) && i >= 0)
- call->unseen = j, call->a[j++] = qsum[i]&3;
- call->n_alleles = j;
+ if (ref_base >= 0) { // for SNPs, find the "unseen" base
+ if (((ref4 < 4 && j < 4) || (ref4 == 4 && j < 5)) && i >= 0)
+ call->unseen = j, call->a[j++] = qsum[i]&3;
+ call->n_alleles = j;
+ } else call->n_alleles = j;
// set the PL array
if (call->n < n) {
call->n = n;
PL[j] = y;
}
}
+// if (ref_base < 0) fprintf(stderr, "%d,%d,%f,%d\n", call->n_alleles, x, sum_min, call->unseen);
call->shift = (int)(sum_min + .499);
}
// combine annotations
b->tid = tid; b->pos = pos; b->qual = 0;
s.s = b->str; s.m = b->m_str; s.l = 0;
kputc('\0', &s);
- kputc("ACGTN"[bc->ori_ref], &s); kputc('\0', &s);
- for (i = 1; i < 5; ++i) {
- if (bc->a[i] < 0) break;
- if (i > 1) kputc(',', &s);
- kputc(bc->unseen == i? 'X' : "ACGT"[bc->a[i]], &s);
+ if (bc->ori_ref < 0) {
+ kputc('N', &s); kputc('\0', &s);
+ kputs("<INDEL>", &s); kputc('\0', &s);
+ } else {
+ kputc("ACGTN"[bc->ori_ref], &s); kputc('\0', &s);
+ for (i = 1; i < 5; ++i) {
+ if (bc->a[i] < 0) break;
+ if (i > 1) kputc(',', &s);
+ kputc(bc->unseen == i? 'X' : "ACGT"[bc->a[i]], &s);
+ }
+ kputc('\0', &s);
}
kputc('\0', &s);
- kputc('\0', &s);
// INFO
kputs("I16=", &s);
for (i = 0; i < 16; ++i) {
bcf_callaux_t *bcf_call_init(double theta, int min_baseQ);
void bcf_call_destroy(bcf_callaux_t *bca);
- int bcf_call_glfgen(int _n, const bam_pileup1_t *pl, int ref_base /*4-bit*/, bcf_callaux_t *bca, bcf_callret1_t *r);
+ int bcf_call_glfgen(int _n, const bam_pileup1_t *pl, int ref_base, bcf_callaux_t *bca, bcf_callret1_t *r);
int bcf_call_combine(int n, const bcf_callret1_t *calls, int ref_base /*4-bit*/, bcf_call_t *call);
int bcf_call2bcf(int tid, int pos, bcf_call_t *bc, bcf1_t *b, bcf_callret1_t *bcr, int is_SP);
int bcf_call_gap_prep(int n, int *n_plp, bam_pileup1_t **plp, int pos, bcf_callaux_t *bca, const char *ref);
#include "ksort.h"
#include "kaln.h"
-#define INDEL_DEBUG
-
#define MINUS_CONST 0x10000000
#define INDEL_WINDOW_SIZE 50
-#define INDEL_BAD_SCORE 10000
static int tpos2qpos(const bam1_core_t *c, const uint32_t *cigar, int32_t tpos, int is_left, int32_t *_tpos)
{
*_tpos = x;
return last_y;
}
+// FIXME: check if the inserted sequence is consistent with the homopolymer run
// l is the relative gap length and l_run is the length of the homopolymer on the reference
static inline int est_seqQ(const bcf_callaux_t *bca, int l, int l_run)
{
sc = ka_global_score((uint8_t*)ref2 + tbeg - left, tend - tbeg + abs(types[t]),
(uint8_t*)query + qbeg, qend - qbeg, &ap);
score[K*n_types + t] = -sc;
-///*
+/*
for (l = 0; l < tend - tbeg + abs(types[t]); ++l)
fputc("ACGTN"[(int)ref2[tbeg-left+l]], stderr);
fputc('\n', stderr);
for (l = 0; l < qend - qbeg; ++l) fputc("ACGTN"[(int)query[qbeg + l]], stderr);
fputc('\n', stderr);
fprintf(stderr, "pos=%d type=%d read=%d:%d name=%s score=%d\n", pos, types[t], s, i, bam1_qname(p->b), sc);
-//*/
+*/
}
}
}
tmp = est_seqQ(bca, types[sc[0]&0x3f], l_run);
}
if (indelQ > tmp) indelQ = tmp;
- if (indelQ > p->b->core.qual) indelQ = p->b->core.qual;
- if (indelQ > bca->capQ) indelQ = bca->capQ;
p->aux = (sc[0]&0x3f)<<8 | indelQ;
sumq[sc[0]&0x3f] += indelQ;
- fprintf(stderr, "pos=%d read=%d:%d name=%s call=%d q=%d,%d\n", pos, s, i, bam1_qname(p->b),
- types[sc[0]&0x3f], indelQ, tmp);
+// 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);
}
}
// determine bca->indel_types[]
for (t = 0; t < 4; ++t) bca->indel_types[t] = B2B_INDEL_NULL;
for (t = 0; t < 4 && t < n_types; ++t)
bca->indel_types[t] = types[sumq[t]&0x3f];
+ // update p->aux
+ for (s = K = 0; s < n; ++s) {
+ for (i = 0; i < n_plp[s]; ++i, ++K) {
+ bam_pileup1_t *p = plp[s] + i;
+ int x = types[p->aux>>8&0x3f];
+ for (j = 0; j < 4; ++j)
+ if (x == bca->indel_types[j]) break;
+ p->aux = j<<8 | (j == 4? 0 : (p->aux&0xff));
+// fprintf(stderr, "pos=%d read=%d:%d name=%s call=%d q=%d\n", pos, s, i, bam1_qname(p->b), p->aux>>8&63, p->aux&0xff);
+ }
+ }
}
// FIXME: to set the inserted sequence
free(score);
bcf_call2bcf(tid, pos, &bc, b, (conf->flag&(MPLP_FMT_DP|MPLP_FMT_SP))? bcr : 0, (conf->flag&MPLP_FMT_SP));
bcf_write(bp, bh, b);
bcf_destroy(b);
-// bcf_call_gap_prep(gplp.n, gplp.n_plp, gplp.plp, pos, bca, ref);
+ // call indels
+ if (bcf_call_gap_prep(gplp.n, gplp.n_plp, gplp.plp, pos, bca, ref) >= 0) {
+ for (i = 0; i < gplp.n; ++i)
+ bcf_call_glfgen(gplp.n_plp[i], gplp.plp[i], -1, bca, bcr + i);
+ bcf_call_combine(gplp.n, bcr, -1, &bc);
+ b = calloc(1, sizeof(bcf1_t));
+ bcf_call2bcf(tid, pos, &bc, b, (conf->flag&(MPLP_FMT_DP|MPLP_FMT_SP))? bcr : 0,
+ (conf->flag&MPLP_FMT_SP));
+ bcf_write(bp, bh, b);
+ bcf_destroy(b);
+ }
} else {
printf("%s\t%d\t%c", h->target_name[tid], pos + 1, (ref && pos < ref_len)? ref[pos] : 'N');
for (i = 0; i < n; ++i) {
ka_param_t ka_param_blast = { 5, 2, 5, 2, aln_sm_blast, 5, 50 };
ka_param_t ka_param_aa2aa = { 10, 2, 10, 2, aln_sm_blosum62, 22, 50 };
-ka_param2_t ka_param2_qual = { 35, 5, 35, 5, 35, 5, 0, 0, aln_sm_qual, 5, 50 };
+ka_param2_t ka_param2_qual = { 37, 11, 37, 11, 37, 11, 0, 0, aln_sm_qual, 5, 50 };
static uint32_t *ka_path2cigar32(const path_t *path, int path_len, int *n_cigar)
{
projects / samtools.git / commitdiff