* make tview more friendly

 * a temporary remedy for an issue in indel calling

diff --git a/bam_maqcns.c b/bam_maqcns.c
index 9a25380f36b77411ecc34336b5a92fe2d2f734d4..cad63d772268db289467c6a4fd625e50b7ac3889 100644 (file)
--- 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 fa5489d13caad7c5c20489091302bfb5b2a5799d..6cc5355fc2db2b1b19e6a3d2d330306a240fa15a 100644 (file)
--- 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 84838dd98e1d56c7270ae78225a8023e9c96ca57..665304f9019621a8734c0651db57a27a26c49fb2 100644 (file)
--- 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 4c121e7cbfa012ac0c3fa61d574f03c893e6d271..7b326fc40e7bf276b7adc3848c0667d5628a68ad 100644 (file)
--- 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;