From 917ed8ca3a137ff2c2db887c0c72f3cce09f7df4 Mon Sep 17 00:00:00 2001 From: Heng Li Date: Mon, 8 Nov 2010 15:54:46 +0000 Subject: [PATCH] indel calling is apparently working, but more information needs to be collected --- bam2bcf.c | 62 +++++++++++++++++++++++++++++++------------------ bam2bcf.h | 2 +- bam2bcf_indel.c | 24 ++++++++++++------- bam_plcmd.c | 12 +++++++++- kaln.c | 2 +- 5 files changed, 68 insertions(+), 34 deletions(-) diff --git a/bam2bcf.c b/bam2bcf.c index 89d6f90..de6a3ed 100644 --- a/bam2bcf.c +++ b/bam2bcf.c @@ -20,7 +20,7 @@ bcf_callaux_t *bcf_call_init(double theta, int min_baseQ) 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; @@ -32,14 +32,17 @@ void bcf_call_destroy(bcf_callaux_t *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; @@ -52,19 +55,24 @@ int bcf_call_glfgen(int _n, const bam_pileup1_t *pl, int ref_base /*4-bit*/, bcf 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; @@ -86,8 +94,10 @@ int bcf_call_combine(int n, const bcf_callret1_t *calls, int ref_base /*4-bit*/, { 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) @@ -108,9 +118,11 @@ int bcf_call_combine(int n, const bcf_callret1_t *calls, int ref_base /*4-bit*/, 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; @@ -138,6 +150,7 @@ int bcf_call_combine(int n, const bcf_callret1_t *calls, int ref_base /*4-bit*/, 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 @@ -158,14 +171,19 @@ int bcf_call2bcf(int tid, int pos, bcf_call_t *bc, bcf1_t *b, bcf_callret1_t *bc 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("", &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) { diff --git a/bam2bcf.h b/bam2bcf.h index 94207d0..e39ab7b 100644 --- a/bam2bcf.h +++ b/bam2bcf.h @@ -36,7 +36,7 @@ extern "C" { 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); diff --git a/bam2bcf_indel.c b/bam2bcf_indel.c index 7ae6be0..1e8a204 100644 --- a/bam2bcf_indel.c +++ b/bam2bcf_indel.c @@ -4,11 +4,8 @@ #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) { @@ -37,6 +34,7 @@ static int tpos2qpos(const bam1_core_t *c, const uint32_t *cigar, 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) { @@ -189,14 +187,14 @@ int bcf_call_gap_prep(int n, int *n_plp, bam_pileup1_t **plp, int pos, bcf_calla 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); -//*/ +*/ } } } @@ -230,12 +228,9 @@ int bcf_call_gap_prep(int n, int *n_plp, bam_pileup1_t **plp, int pos, bcf_calla 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[] @@ -254,6 +249,17 @@ int bcf_call_gap_prep(int n, int *n_plp, bam_pileup1_t **plp, int pos, bcf_calla 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); diff --git a/bam_plcmd.c b/bam_plcmd.c index 1b15aae..b56f879 100644 --- a/bam_plcmd.c +++ b/bam_plcmd.c @@ -729,7 +729,17 @@ static int mpileup(mplp_conf_t *conf, int n, char **fn) 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) { diff --git a/kaln.c b/kaln.c index 01efe2e..9c0bbaa 100644 --- a/kaln.c +++ b/kaln.c @@ -84,7 +84,7 @@ int aln_sm_qual[] = { 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) { -- 2.39.2