From 4c8c9dfc1e3b3b066a62a703fd3ba04db6ad5a45 Mon Sep 17 00:00:00 2001 From: Heng Li Date: Thu, 18 Nov 2010 02:41:12 +0000 Subject: [PATCH] * samtools-0.1.10 (r831) * initial MNP caller --- bam2bcf.c | 33 +++++-- bam2bcf.h | 9 +- bam2bcf_indel.c | 255 ++++++++++++++++++++++++++++++++++++++++++++++++ bam_plcmd.c | 35 +++++-- bamtk.c | 2 +- 5 files changed, 313 insertions(+), 21 deletions(-) diff --git a/bam2bcf.c b/bam2bcf.c index 088635c..fe447e6 100644 --- a/bam2bcf.c +++ b/bam2bcf.c @@ -22,6 +22,7 @@ bcf_callaux_t *bcf_call_init(double theta, int min_baseQ) bca->capQ = 60; bca->openQ = 40; bca->extQ = 20; bca->tandemQ = 100; bca->min_baseQ = min_baseQ; + bca->last_mnp_pos = -100; bca->e = errmod_init(1. - theta); return bca; } @@ -32,8 +33,10 @@ void bcf_call_destroy(bcf_callaux_t *bca) errmod_destroy(bca->e); free(bca->bases); free(bca->inscns); free(bca); } -/* ref_base is the 4-bit representation of the reference base. It is - * negative if we are looking at an indel. */ +/* Compute genotype likelihood by combining likelihood of each + * read. ref_base is the 4-bit representation of the reference base. It + * is negative if we are looking at an indel or an MNP. Indel and MNP + * are indistinguishable in this function. */ 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, is_indel, ori_depth = 0; @@ -92,7 +95,9 @@ int bcf_call_glfgen(int _n, const bam_pileup1_t *pl, int ref_base, bcf_callaux_t errmod_cal(bca->e, n, 5, bca->bases, r->p); return r->depth; } - +/* Combine individual calls. ref_base is a 4-bit A/C/G/T for a SNP, or + * B2B_REF_MNP, or B2B_REF_INDEL. Indel and MNP are indistinguishable in + * this function. */ int bcf_call_combine(int n, const bcf_callret1_t *calls, int ref_base /*4-bit*/, bcf_call_t *call) { int ref4, i, j, qsum[4]; @@ -100,7 +105,7 @@ int bcf_call_combine(int n, const bcf_callret1_t *calls, int ref_base /*4-bit*/, 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; + } else call->ori_ref = ref_base, ref4 = 0; // ref_base < 0 // calculate qsum memset(qsum, 0, 4 * sizeof(int)); for (i = 0; i < n; ++i) @@ -125,7 +130,7 @@ int bcf_call_combine(int n, const bcf_callret1_t *calls, int ref_base /*4-bit*/, 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 { + } else { // for indels and MNPs call->n_alleles = j; if (call->n_alleles == 1) return -1; // no reliable supporting read. stop doing anything } @@ -168,7 +173,8 @@ int bcf_call_combine(int n, const bcf_callret1_t *calls, int ref_base /*4-bit*/, } return 0; } - +/* Fill a bcf1_t record. The type of the record (SNP, MNP or INDEL) is + * determined by bca->ori_ref. */ int bcf_call2bcf(int tid, int pos, bcf_call_t *bc, bcf1_t *b, bcf_callret1_t *bcr, int is_SP, const bcf_callaux_t *bca, const char *ref) { @@ -179,7 +185,7 @@ 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); - if (bc->ori_ref < 0) { // an indel + if (bc->ori_ref == B2B_REF_INDEL) { // an indel // write REF kputc(ref[pos], &s); for (j = 0; j < bca->indelreg; ++j) kputc(ref[pos+1+j], &s); @@ -202,6 +208,16 @@ int bcf_call2bcf(int tid, int pos, bcf_call_t *bc, bcf1_t *b, bcf_callret1_t *bc } } kputc('\0', &s); + } else if (bc->ori_ref == B2B_REF_MNP) { + for (j = 0; j < bca->indelreg; ++j) kputc(ref[pos+1+j], &s); + kputc('\0', &s); + for (i = 1; i < 4; ++i) { + if (bc->a[i] < 0) break; + if (i > 1) kputc(',', &s); + for (j = 0; j < bca->indelreg; ++j) + kputc("ACGT"[bca->indel_types[i]>>2*i&3], &s); + } + kputc('\0', &s); } else { // a SNP kputc("ACGTN"[bc->ori_ref], &s); kputc('\0', &s); for (i = 1; i < 5; ++i) { @@ -213,7 +229,8 @@ int bcf_call2bcf(int tid, int pos, bcf_call_t *bc, bcf1_t *b, bcf_callret1_t *bc } kputc('\0', &s); // INFO - if (bc->ori_ref < 0) kputs("INDEL;", &s); + if (bc->ori_ref == B2B_REF_INDEL) kputs("INDEL;", &s); + else if (bc->ori_ref == B2B_REF_MNP) kputs("MNP;", &s); kputs("DP=", &s); kputw(bc->ori_depth, &s); kputs(";I16=", &s); for (i = 0; i < 16; ++i) { if (i) kputc(',', &s); diff --git a/bam2bcf.h b/bam2bcf.h index 26b022c..0de76f9 100644 --- a/bam2bcf.h +++ b/bam2bcf.h @@ -6,6 +6,10 @@ #include "bcftools/bcf.h" #define B2B_INDEL_NULL 10000 +#define B2B_MAX_MNP 4 // cannot be larger than 4!!! +#define B2B_MNP_WIN 10 +#define B2B_REF_INDEL (-1) +#define B2B_REF_MNP (-2) typedef struct __bcf_callaux_t { int capQ, min_baseQ; @@ -13,7 +17,7 @@ typedef struct __bcf_callaux_t { // for internal uses int max_bases; int indel_types[4]; - int maxins, indelreg; + int maxins, indelreg, last_mnp_pos; char *inscns; uint16_t *bases; errmod_t *e; @@ -28,7 +32,7 @@ typedef struct { typedef struct { int a[5]; // alleles: ref, alt, alt2, alt3 - int n, n_alleles, shift, ori_ref, unseen; + int n, n_alleles, shift, ori_ref, unseen; // ori_ref can be B2B_REF_INDEL/B2B_REF_MNP int anno[16], depth, ori_depth; uint8_t *PL; } bcf_call_t; @@ -45,6 +49,7 @@ extern "C" { const bcf_callaux_t *bca, const char *ref); int bcf_call_gap_prep(int n, int *n_plp, bam_pileup1_t **plp, int pos, bcf_callaux_t *bca, const char *ref, const void *rghash); + int bcf_call_mnp_prep(int n, int *n_plp, bam_pileup1_t **plp, int pos, bcf_callaux_t *bca, const char *ref); #ifdef __cplusplus } diff --git a/bam2bcf_indel.c b/bam2bcf_indel.c index 819c1e7..8b80e3f 100644 --- a/bam2bcf_indel.c +++ b/bam2bcf_indel.c @@ -408,3 +408,258 @@ int bcf_call_gap_prep(int n, int *n_plp, bam_pileup1_t **plp, int pos, bcf_calla free(types); free(inscns); return n_alt > 0? 0 : -1; } + +#define END_SKIP 4 // must be larger than B2B_MAX_MNP +#define MIN_MNP_FLANK_BAQ 20 + +// the following routine is modified from bcf_call_gap_prep() +int bcf_call_mnp_prep(int n, int *n_plp, bam_pileup1_t **plp, int pos, bcf_callaux_t *bca, const char *ref) +{ + extern void ks_introsort_uint32_t(int, uint32_t*); + int i, s, j, k, t, n_types, *types, max_rd_len, left, right, *score1; + int N, K, ref_seq, ref_type, n_alt; + char *ref2, *query; + uint8_t *_seqQ; + if (ref == 0 || bca == 0) return -1; + if (pos < bca->last_mnp_pos + B2B_MNP_WIN) return -2; // to avoid calling a TNP multiple times + if (pos < END_SKIP || ref[pos] == 0 || ref[pos+1] == 0) return -2; // end of the reference + { // determine if there is an MNP + int r[2]; + uint8_t *tt; + r[0] = bam_nt16_table[(int)ref[pos]]; + r[1] = bam_nt16_table[(int)ref[pos+1]]; + for (s = 0; s < n; ++s) { + for (i = 0; i < n_plp[s]; ++i) { + bam_pileup1_t *p = plp[s] + i; + int left, rght; + if (p->qpos < END_SKIP || p->qpos >= p->b->core.l_qseq - END_SKIP) continue; + tt = bam1_seq(p->b); + if (bam1_seqi(tt, p->qpos) == r[0] || bam1_seqi(tt, p->qpos+1) == r[1]) continue; // no MNP + tt = bam1_qual(p->b); + for (left = 0, k = p->qpos - END_SKIP; k <= p->qpos; ++k) + left = left > tt[k]? left : tt[k]; + for (rght = 0, k = p->qpos; k < p->b->core.l_qseq - END_SKIP; ++k) + rght = rght > tt[k]? left : tt[k]; + if (left >= MIN_MNP_FLANK_BAQ && rght >= MIN_MNP_FLANK_BAQ) break; // bracketed by good bases + } + if (i != n_plp[s]) break; + } + if (s == n) return -1; // there is no MNP at this position. + } + for (s = N = 0; s < n; ++s) N += n_plp[s]; // N is the total number of reads + { // construct the MNP consensus + uint8_t *tt; + uint32_t *aux, x, *cnt; + int m; + m = max_rd_len = 0; + aux = calloc(N + 1, 4); + for (i = 0, x = 0; i < B2B_MAX_MNP; ++i) + x |= bam_nt16_nt4_table[bam_nt16_table[(int)ref[pos+i]]] << 2*i; + ref_seq = aux[m++] = x; + bca->indelreg = 0; + for (s = 0; s < n; ++s) { + for (i = 0; i < n_plp[s]; ++i) { + bam_pileup1_t *p = plp[s] + i; + int stop; + j = bam_cigar2qlen(&p->b->core, bam1_cigar(p->b)); + if (j > max_rd_len) max_rd_len = j; + if (p->qpos < END_SKIP || p->qpos >= p->b->core.l_qseq - END_SKIP) continue; + tt = bam1_seq(p->b); + for (k = j = stop = 0, x = 0; k < B2B_MAX_MNP; ++k) { + int c = bam_nt16_nt4_table[bam1_seqi(tt, p->qpos + k)]; + if (c > 3) break; + if (c != (ref_seq>>k*2&3) && !stop) ++j; + else stop = 1; + x |= c << 2*k; + } +// fprintf(stderr, "* %d, %d\n", x, j); + if (k == B2B_MAX_MNP && j >= 2) aux[m++] = x; + } + } + ks_introsort(uint32_t, m, aux); + // squeeze out indentical types + for (i = 1, n_types = 1; i < m; ++i) + if (aux[i] != aux[i-1]) ++n_types; + if (n_types == 1) { // no indels + free(aux); return -1; + } + // count reads for each type + cnt = alloca(n_types * 4); + cnt[0] = 1<<8 | aux[0]; + for (i = 1, t = 0; i < m; ++i) { + if (aux[i] != aux[i-1]) { + ++t; + cnt[t] = 1<<8 | aux[i]; + } else cnt[t] += 1<<8; + } + free(aux); + // collect types + ks_introsort(uint32_t, n_types, cnt); + for (t = 0; t < n_types; ++t) + if ((cnt[t]&0xff) == ref_seq) break; + if (t == n_types - 1) --t; // we get at least two types + types = (int*)calloc(n_types, sizeof(int)); + for (i = 0; t < n_types; ++i, ++t) types[i] = cnt[t]&0xff; + n_types = i; + // calculate MNP length + for (i = B2B_MAX_MNP - 1, k = 0; i >= 0; --i) { + int c = types[0] >> 2*i & 3; + for (t = 1; t < n_types; ++t) + if ((types[t] >> 2*i & 3) != c) break; + if (t == n_types) ++k; + else break; + } + bca->indelreg = B2B_MAX_MNP - k; + x = (1<<2*bca->indelreg) - 1; + for (t = 0; t < n_types; ++t) types[t] &= x; + ref_seq &= x; + for (t = 0; t < n_types; ++t) + if (types[t] == ref_seq) break; + ref_type = t; +// fprintf(stderr, "%d, %d, [%x,%x] %d\n", pos, n_types, types[0], types[1], bca->indelreg); + } + { // calculate left and right boundary + left = pos > INDEL_WINDOW_SIZE? pos - INDEL_WINDOW_SIZE : 0; + right = pos + INDEL_WINDOW_SIZE; + // in case the alignments stand out the reference + for (i = pos; i < right; ++i) + if (ref[i] == 0) break; + right = i; + } + { // find the highest base quality + _seqQ = calloc(N, 1); + for (s = K = 0; s < n; ++s) { + for (i = 0; i < n_plp[s]; ++i, ++K) { + bam_pileup1_t *p = plp[s] + i; + uint8_t *qq, *bq; + int max = 0; + qq = bam1_qual(p->b); + bq = (uint8_t*)bam_aux_get(p->b, "ZQ"); + if (bq) ++bq; + for (j = p->qpos; j < p->b->core.l_qseq && j < p->qpos + bca->indelreg; ++j) { + int q = bq? qq[j] + (bq[j] - 64) : qq[j]; + max = max > q? max : q; + } + _seqQ[K] = max; + } + } + } + // compute the likelihood given each type of indel for each read + ref2 = calloc(right - left + 2, 1); + query = calloc(right - left + max_rd_len + 2, 1); + score1 = calloc(N * n_types, sizeof(int)); + for (t = 0; t < n_types; ++t) { + int l; + kpa_par_t apf1 = { 1e-4, 1e-2, 10 }; + apf1.bw = abs(types[t]) + 3; + // write ref2 + for (k = 0, j = left; j < pos; ++j) + ref2[k++] = bam_nt16_nt4_table[bam_nt16_table[(int)ref[j]]]; + for (i = 0; i < bca->indelreg; ++i, ++j) + ref2[k++] = types[t]>>2*i&3; + for (; j < right && ref[j]; ++j) + ref2[k++] = bam_nt16_nt4_table[bam_nt16_table[(int)ref[j]]]; + if (j < right) right = j; + // align each read to ref2 + for (s = K = 0; s < n; ++s) { + for (i = 0; i < n_plp[s]; ++i, ++K) { + bam_pileup1_t *p = plp[s] + i; + int qbeg, qend, tbeg, tend, sc; + uint8_t *seq = bam1_seq(p->b); + // determine the start and end of sequences for alignment + qbeg = tpos2qpos(&p->b->core, bam1_cigar(p->b), left, 0, &tbeg); + qend = tpos2qpos(&p->b->core, bam1_cigar(p->b), right, 1, &tend); + // write the query sequence + for (l = qbeg; l < qend; ++l) + query[l - qbeg] = bam_nt16_nt4_table[bam1_seqi(seq, l)]; + { // do realignment; this is the bottleneck + const uint8_t *qual = bam1_qual(p->b), *bq; + uint8_t *qq; + qq = calloc(qend - qbeg, 1); + bq = (uint8_t*)bam_aux_get(p->b, "ZQ"); + if (bq) ++bq; // skip type + for (l = qbeg; l < qend; ++l) { + qq[l - qbeg] = bq? qual[l] + (bq[l] - 64) : qual[l]; + if (qq[l - qbeg] > 30) qq[l - qbeg] = 30; + if (qq[l - qbeg] < 7) qq[l - qbeg] = 7; + } + sc = kpa_glocal((uint8_t*)ref2 + tbeg - left, tend - tbeg, + (uint8_t*)query, qend - qbeg, qq, &apf1, 0, 0); + l = (int)(100. * sc / (qend - qbeg) + .499); // used for adjusting indelQ below + if (l > 255) l = 255; + score1[K*n_types + t] = sc<<8 | l; + free(qq); + } +/* + for (l = 0; l < tend - tbeg; ++l) + fputc("ACGTN"[(int)ref2[tbeg-left+l]], stderr); + fputc('\n', stderr); + for (l = 0; l < qend - qbeg; ++l) fputc("ACGTN"[(int)query[l]], stderr); + fputc('\n', stderr); + fprintf(stderr, "pos=%d type=%d read=%d:%d name=%s qbeg=%d tbeg=%d score=%d\n", pos, types[t], s, i, bam1_qname(p->b), qbeg, tbeg, sc); +*/ + } + } + } + free(ref2); free(query); + { // compute indelQ + int *sc, tmp, *sumq; + sc = alloca(n_types * sizeof(int)); + sumq = alloca(n_types * sizeof(int)); + memset(sumq, 0, sizeof(int) * n_types); + for (s = K = 0; s < n; ++s) { + for (i = 0; i < n_plp[s]; ++i, ++K) { + bam_pileup1_t *p = plp[s] + i; + int *sct = &score1[K*n_types], seqQ, indelQ; + for (t = 0; t < n_types; ++t) sc[t] = sct[t]<<6 | t; + for (t = 1; t < n_types; ++t) // insertion sort + for (j = t; j > 0 && sc[j] < sc[j-1]; --j) + tmp = sc[j], sc[j] = sc[j-1], sc[j-1] = tmp; + if ((sc[0]&0x3f) == ref_type) { + indelQ = (sc[1]>>14) - (sc[0]>>14); + } else { + for (t = 0; t < n_types; ++t) // look for the reference type + if ((sc[t]&0x3f) == ref_type) break; + indelQ = (sc[t]>>14) - (sc[0]>>14); + } + seqQ = _seqQ[K]; + tmp = sc[0]>>6 & 0xff; + indelQ = tmp > 111? 0 : (int)((1. - tmp/111.) * indelQ + .499); // reduce indelQ + p->aux = (sc[0]&0x3f)<<16 | seqQ<<8 | indelQ; + sumq[sc[0]&0x3f] += indelQ < seqQ? indelQ : seqQ; + } + } + // determine bca->indel_types[] and bca->inscns + for (t = 0; t < n_types; ++t) + sumq[t] = sumq[t]<<6 | t; + for (t = 1; t < n_types; ++t) // insertion sort + for (j = t; j > 0 && sumq[j] > sumq[j-1]; --j) + tmp = sumq[j], sumq[j] = sumq[j-1], sumq[j-1] = tmp; + for (t = 0; t < n_types; ++t) // look for the reference type + if ((sumq[t]&0x3f) == ref_type) break; + if (t) { // then move the reference type to the first + tmp = sumq[t]; + for (; t > 0; --t) sumq[t] = sumq[t-1]; + sumq[0] = tmp; + } + 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 = n_alt = 0; s < n; ++s) { + for (i = 0; i < n_plp[s]; ++i) { + bam_pileup1_t *p = plp[s] + i; + int x = types[p->aux>>16&0x3f]; + for (j = 0; j < 4; ++j) + if (x == bca->indel_types[j]) break; + p->aux = j<<16 | (j == 4? 0 : (p->aux&0xffff)); + if ((p->aux>>16&0x3f) > 0) ++n_alt; +// fprintf(stderr, "X pos=%d read=%d:%d name=%s call=%d type=%d q=%d seqQ=%d\n", pos, s, i, bam1_qname(p->b), p->aux>>16&63, bca->indel_types[p->aux>>16&63], p->aux&0xff, p->aux>>8&0xff); + } + } + } + free(score1); + // free + free(types); free(_seqQ); + return n_alt > 0? 0 : -1; +} diff --git a/bam_plcmd.c b/bam_plcmd.c index e3f73aa..7f13eee 100644 --- a/bam_plcmd.c +++ b/bam_plcmd.c @@ -738,16 +738,31 @@ static int mpileup(mplp_conf_t *conf, int n, char **fn) (conf->flag&MPLP_FMT_SP), 0, 0); bcf_write(bp, bh, b); bcf_destroy(b); - // call indels - if (!(conf->flag&MPLP_NO_INDEL) && bcf_call_gap_prep(gplp.n, gplp.n_plp, gplp.plp, pos, bca, ref, rghash) >= 0) { - for (i = 0; i < gplp.n; ++i) - bcf_call_glfgen(gplp.n_plp[i], gplp.plp[i], -1, bca, bcr + i); - if (bcf_call_combine(gplp.n, bcr, -1, &bc) >= 0) { - 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), bca, ref); - bcf_write(bp, bh, b); - bcf_destroy(b); + if (!(conf->flag&MPLP_NO_INDEL)) { + // call MNPs + if (bcf_call_mnp_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], B2B_REF_MNP, bca, bcr + i); + if (bcf_call_combine(gplp.n, bcr, B2B_REF_MNP, &bc) >= 0) { + 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), bca, ref); + bcf_write(bp, bh, b); + bcf_destroy(b); + bca->last_mnp_pos = pos; + } + } + // call indels + if (bcf_call_gap_prep(gplp.n, gplp.n_plp, gplp.plp, pos, bca, ref, rghash) >= 0) { + for (i = 0; i < gplp.n; ++i) + bcf_call_glfgen(gplp.n_plp[i], gplp.plp[i], B2B_REF_INDEL, bca, bcr + i); + if (bcf_call_combine(gplp.n, bcr, B2B_REF_INDEL, &bc) >= 0) { + 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), bca, ref); + bcf_write(bp, bh, b); + bcf_destroy(b); + } } } } else { diff --git a/bamtk.c b/bamtk.c index 88328d4..6d7b832 100644 --- a/bamtk.c +++ b/bamtk.c @@ -9,7 +9,7 @@ #endif #ifndef PACKAGE_VERSION -#define PACKAGE_VERSION "0.1.10 (r829)" +#define PACKAGE_VERSION "0.1.10-1 (r831)" #endif int bam_taf2baf(int argc, char *argv[]); -- 2.39.2