X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;f=bam_tview.c;h=20e6f466dfdf9ab0a2ba444e4181a20115b0632e;hb=9f118264ea012adc21a46d7c03eaad4f9ce7d4d4;hp=e48afa7bd4add26404fecc04d6548e7b3989927b;hpb=253dbf083d8a84721ae6ac63af786c8dcf566ecb;p=samtools.git diff --git a/bam_tview.c b/bam_tview.c index e48afa7..20e6f46 100644 --- a/bam_tview.c +++ b/bam_tview.c @@ -19,9 +19,12 @@ #include #include #include +#include +#include #include "bam.h" #include "faidx.h" -#include "bam_maqcns.h" +#include "bam2bcf.h" +#include "sam_header.h" char bam_aux_getCEi(bam1_t *b, int i); char bam_aux_getCSi(bam1_t *b, int i); @@ -50,14 +53,17 @@ typedef struct { bamFile fp; int curr_tid, left_pos; faidx_t *fai; - bam_maqcns_t *bmc; + bcf_callaux_t *bca; int ccol, last_pos, row_shift, base_for, color_for, is_dot, l_ref, ins, no_skip, show_name; char *ref; + char *sample; //TODO: multiple samples and read groups + void *rg2sm; } tview_t; int tv_pl_func(uint32_t tid, uint32_t pos, int n, const bam_pileup1_t *pl, void *data) { + extern unsigned char bam_nt16_table[256]; tview_t *tv = (tview_t*)data; int i, j, c, rb, attr, max_ins = 0; uint32_t call = 0; @@ -70,11 +76,26 @@ int tv_pl_func(uint32_t tid, uint32_t pos, int n, const bam_pileup1_t *pl, void mvaddch(1, tv->ccol++, c); } if (pos%10 == 0 && tv->mcol - tv->ccol >= 10) mvprintw(0, tv->ccol, "%-d", pos+1); - // print consensus - call = bam_maqcns_call(n, pl, tv->bmc); + { // call consensus + bcf_callret1_t bcr; + int qsum[4], a1, a2, tmp; + double p[3], prior = 30; + bcf_call_glfgen(n, pl, bam_nt16_table[rb], tv->bca, &bcr); + for (i = 0; i < 4; ++i) qsum[i] = bcr.qsum[i]<<2 | i; + for (i = 1; i < 4; ++i) // insertion sort + for (j = i; j > 0 && qsum[j] > qsum[j-1]; --j) + tmp = qsum[j], qsum[j] = qsum[j-1], qsum[j-1] = tmp; + a1 = qsum[0]&3; a2 = qsum[1]&3; + p[0] = bcr.p[a1*5+a1]; p[1] = bcr.p[a1*5+a2] + prior; p[2] = bcr.p[a2*5+a2]; + if ("ACGT"[a1] != toupper(rb)) p[0] += prior + 3; + if ("ACGT"[a2] != toupper(rb)) p[2] += prior + 3; + if (p[0] < p[1] && p[0] < p[2]) call = (1<>28&0xf]; - i = (call>>8&0xff)/10+1; + c = ",ACMGRSVTWYHKDBN"[call>>16&0xf]; + i = (call&0xffff)/10+1; if (i > 4) i = 4; attr |= COLOR_PAIR(i); if (c == toupper(rb)) c = '.'; @@ -97,7 +118,6 @@ int tv_pl_func(uint32_t tid, uint32_t pos, int n, const bam_pileup1_t *pl, void if (!p->is_del) { if (tv->base_for == TV_BASE_COLOR_SPACE && (c = bam_aux_getCSi(p->b, p->qpos))) { - c = bam_aux_getCSi(p->b, p->qpos); // assume that if we found one color, we will be able to get the color error if (tv->is_dot && '-' == bam_aux_getCEi(p->b, p->qpos)) c = bam1_strand(p->b)? ',' : '.'; } else { @@ -179,21 +199,27 @@ int tv_pl_func(uint32_t tid, uint32_t pos, int n, const bam_pileup1_t *pl, void return 0; } -tview_t *tv_init(const char *fn, const char *fn_fa) +tview_t *tv_init(const char *fn, const char *fn_fa, char *samples) { tview_t *tv = (tview_t*)calloc(1, sizeof(tview_t)); tv->is_dot = 1; - tv->idx = bam_index_load(fn); - if (tv->idx == 0) exit(1); tv->fp = bam_open(fn, "r"); bgzf_set_cache_size(tv->fp, 8 * 1024 *1024); assert(tv->fp); tv->header = bam_header_read(tv->fp); + tv->idx = bam_index_load(fn); + if (tv->idx == 0) exit(1); tv->lplbuf = bam_lplbuf_init(tv_pl_func, tv); if (fn_fa) tv->fai = fai_load(fn_fa); - tv->bmc = bam_maqcns_init(); + tv->bca = bcf_call_init(0.83, 13); tv->ins = 1; - bam_maqcns_prepare(tv->bmc); + + if ( samples ) + { + tv->sample = samples; + tv->header->dict = sam_header_parse2(tv->header->text); + tv->rg2sm = sam_header2tbl(tv->header->dict, "RG", "ID", "SM"); + } initscr(); keypad(stdscr, TRUE); @@ -224,7 +250,7 @@ void tv_destroy(tview_t *tv) endwin(); bam_lplbuf_destroy(tv->lplbuf); - bam_maqcns_destroy(tv->bmc); + bcf_call_destroy(tv->bca); bam_index_destroy(tv->idx); if (tv->fai) fai_destroy(tv->fai); free(tv->ref); @@ -236,6 +262,14 @@ void tv_destroy(tview_t *tv) int tv_fetch_func(const bam1_t *b, void *data) { tview_t *tv = (tview_t*)data; + if ( tv->sample ) + { + const uint8_t *rg = bam_aux_get(b, "RG"); + if ( !rg ) return 0; + const char *sm = sam_tbl_get(tv->rg2sm, (const char*)(rg + 1)); + if ( !sm ) return 0; + if ( strcmp(sm,tv->sample) ) return 0; + } if (tv->no_skip) { uint32_t *cigar = bam1_cigar(b); // this is cheating... int i; @@ -288,7 +322,7 @@ static void tv_win_goto(tview_t *tv, int *tid, int *pos) int c = wgetch(tv->wgoto); wrefresh(tv->wgoto); if (c == KEY_BACKSPACE || c == '\010' || c == '\177') { - --l; + if(l > 0) --l; } else if (c == KEY_ENTER || c == '\012' || c == '\015') { int _tid = -1, _beg, _end; if (str[0] == '=') { @@ -400,14 +434,40 @@ end_loop: return; } +void error(const char *format, ...) +{ + if ( !format ) + { + fprintf(stderr, "\n"); + fprintf(stderr, "Usage: bamtk tview [options] [ref.fasta]\n"); + fprintf(stderr, "Options:\n"); + fprintf(stderr, " -s STR display only reads from this sample\n"); + fprintf(stderr, "\n\n"); + } + else + { + va_list ap; + va_start(ap, format); + vfprintf(stderr, format, ap); + va_end(ap); + } + exit(-1); +} + + int bam_tview_main(int argc, char *argv[]) { tview_t *tv; - if (argc == 1) { - fprintf(stderr, "Usage: bamtk tview [ref.fasta]\n"); - return 1; - } - tv = tv_init(argv[1], (argc == 2)? 0 : argv[2]); + char *samples=NULL; + int c; + while ((c = getopt(argc, argv, "s:")) >= 0) { + switch (c) { + case 's': samples=optarg; break; + default: error(NULL); + } + } + if (argc==optind) error(NULL); + tv = tv_init(argv[optind], (optind+1>=argc)? 0 : argv[optind+1], samples); tv_draw_aln(tv, 0, 0); tv_loop(tv); tv_destroy(tv);