From 961c3fe158b5a9777bc8c061425ca3cb15d8687c Mon Sep 17 00:00:00 2001 From: Heng Li Date: Fri, 9 Jul 2010 02:36:12 +0000 Subject: [PATCH] * make tview more friendly * a temporary remedy for an issue in indel calling --- bam_maqcns.c | 7 +++++-- bam_maqcns.h | 5 +++-- bam_plcmd.c | 2 +- bam_tview.c | 19 ++++++++++++++----- 4 files changed, 23 insertions(+), 10 deletions(-) diff --git a/bam_maqcns.c b/bam_maqcns.c index 9a25380..cad63d7 100644 --- a/bam_maqcns.c +++ b/bam_maqcns.c @@ -310,6 +310,7 @@ bam_maqindel_opt_t *bam_maqindel_opt_init() bam_maqindel_opt_t *mi = (bam_maqindel_opt_t*)calloc(1, sizeof(bam_maqindel_opt_t)); mi->q_indel = 40; mi->r_indel = 0.00015; + mi->r_snp = 0.001; // mi->mm_penalty = 3; mi->indel_err = 4; @@ -406,7 +407,8 @@ bam_maqindel_ret_t *bam_maqindel(int n, int pos, const bam_maqindel_opt_t *mi, c } { // the core part char *ref2, *rs, *inscns = 0; - int k, l, *score, *pscore, max_ins = types[n_types-1]; + int qr_snp, k, l, *score, *pscore, max_ins = types[n_types-1]; + qr_snp = (int)(-4.343 * log(mi->r_snp) + .499); if (max_ins > 0) { // get the consensus of inserted sequences int *inscns_aux = (int*)calloc(4 * n_types * max_ins, sizeof(int)); // count occurrences @@ -488,7 +490,8 @@ bam_maqindel_ret_t *bam_maqindel(int n, int pos, const bam_maqindel_opt_t *mi, c if (op == BAM_CMATCH) { int k; for (k = 0; k < len; ++k) - if (ref2[x+k] != rs[y+k] && ref2[x+k] < 4) ps += bam1_qual(p->b)[y+k]; + if (ref2[x+k] != rs[y+k] && ref2[x+k] < 4) + ps += bam1_qual(p->b)[y+k] < qr_snp? bam1_qual(p->b)[y+k] : qr_snp; x += len; y += len; } else if (op == BAM_CINS || op == BAM_CSOFT_CLIP) { if (op == BAM_CINS && l > 0 && l < n_acigar - 1) ps += mi->q_indel * len; diff --git a/bam_maqcns.h b/bam_maqcns.h index fa5489d..6cc5355 100644 --- a/bam_maqcns.h +++ b/bam_maqcns.h @@ -16,8 +16,9 @@ typedef struct { } bam_maqcns_t; typedef struct { - int q_indel; - float r_indel; + int q_indel; // indel sequencing error, phred scaled + float r_indel; // indel prior + float r_snp; // snp prior // hidden parameters, unchangeable from command line int mm_penalty, indel_err, ambi_thres; } bam_maqindel_opt_t; diff --git a/bam_plcmd.c b/bam_plcmd.c index 84838dd..665304f 100644 --- a/bam_plcmd.c +++ b/bam_plcmd.c @@ -352,7 +352,7 @@ int bam_pileup(int argc, char *argv[]) case 'f': fn_fa = strdup(optarg); break; case 'T': d->c->theta = atof(optarg); break; case 'N': d->c->n_hap = atoi(optarg); break; - case 'r': d->c->het_rate = atof(optarg); break; + case 'r': d->c->het_rate = atof(optarg); d->ido->r_snp = d->c->het_rate; break; case 'M': d->c->cap_mapQ = atoi(optarg); break; case 'd': d->max_depth = atoi(optarg); break; case 'c': d->format |= BAM_PLF_CNS; break; diff --git a/bam_tview.c b/bam_tview.c index 4c121e7..7b326fc 100644 --- a/bam_tview.c +++ b/bam_tview.c @@ -280,7 +280,7 @@ int tv_draw_aln(tview_t *tv, int tid, int pos) static void tv_win_goto(tview_t *tv, int *tid, int *pos) { - char str[256]; + char str[256], *p; int i, l = 0; wborder(tv->wgoto, '|', '|', '-', '-', '+', '+', '+', '+'); mvwprintw(tv->wgoto, 1, 2, "Goto: "); @@ -291,10 +291,18 @@ static void tv_win_goto(tview_t *tv, int *tid, int *pos) --l; } else if (c == KEY_ENTER || c == '\012' || c == '\015') { int _tid = -1, _beg, _end; - bam_parse_region(tv->header, str, &_tid, &_beg, &_end); - if (_tid >= 0) { - *tid = _tid; *pos = _beg; - return; + if (str[0] == '=') { + _beg = strtol(str+1, &p, 10); + if (_beg > 0) { + *pos = _beg; + return; + } + } else { + bam_parse_region(tv->header, str, &_tid, &_beg, &_end); + if (_tid >= 0) { + *tid = _tid; *pos = _beg; + return; + } } } else if (isgraph(c)) { if (l < TV_MAX_GOTO) str[l++] = c; @@ -351,6 +359,7 @@ void tv_loop(tview_t *tv) case '?': tv_win_help(tv); break; case '\033': case 'q': goto end_loop; + case '/': case 'g': tv_win_goto(tv, &tid, &pos); break; case 'm': tv->color_for = TV_COLOR_MAPQ; break; case 'b': tv->color_for = TV_COLOR_BASEQ; break; -- 2.39.2