From 5c4ef838b286eeb6fffc7c41ffbb776281963c8e Mon Sep 17 00:00:00 2001
From: Heng Li <lh3@live.co.uk>
Date: Fri, 20 Aug 2010 19:46:40 +0000
Subject: [PATCH] added revised MAQ error model
---
Makefile | 2 +-
bam_maqcns.c | 43 +++++++++++------
bam_maqcns.h | 6 ++-
bam_plcmd.c | 6 +--
errmod.c | 130 +++++++++++++++++++++++++++++++++++++++++++++++++++
errmod.h | 17 +++++++
ksort.h | 10 ++++
7 files changed, 194 insertions(+), 20 deletions(-)
create mode 100644 errmod.c
create mode 100644 errmod.h
diff --git a/Makefile b/Makefile
index 6a04c1b..c14d663 100644
--- a/Makefile
+++ b/Makefile
@@ -7,7 +7,7 @@ LOBJS= bgzf.o kstring.o bam_aux.o bam.o bam_import.o sam.o bam_index.o \
$(KNETFILE_O) bam_sort.o sam_header.o bam_reheader.o
AOBJS= bam_tview.o bam_maqcns.o bam_plcmd.o sam_view.o \
bam_rmdup.o bam_rmdupse.o bam_mate.o bam_stat.o bam_color.o \
- bamtk.o kaln.o bam_mcns.o bam2bcf.o
+ bamtk.o kaln.o bam_mcns.o bam2bcf.o errmod.o
PROG= samtools
INCLUDES= -I.
SUBDIRS= . bcftools misc
diff --git a/bam_maqcns.c b/bam_maqcns.c
index cad63d7..9caa5d7 100644
--- a/bam_maqcns.c
+++ b/bam_maqcns.c
@@ -3,6 +3,7 @@
#include "bam.h"
#include "bam_maqcns.h"
#include "ksort.h"
+#include "errmod.h"
#include "kaln.h"
KSORT_INIT_GENERIC(uint32_t)
@@ -12,12 +13,13 @@ KSORT_INIT_GENERIC(uint32_t)
typedef struct __bmc_aux_t {
int max;
uint32_t *info;
+ uint16_t *info16;
+ errmod_t *em;
} bmc_aux_t;
typedef struct {
float esum[4], fsum[4];
uint32_t c[4];
- uint32_t rms_mapQ;
} glf_call_aux_t;
char bam_nt16_nt4_table[] = { 4, 0, 1, 4, 2, 4, 4, 4, 3, 4, 4, 4, 4, 4, 4, 4 };
@@ -41,7 +43,7 @@ static void cal_het(bam_maqcns_t *aa)
for (n1 = 0; n1 < 256; ++n1) {
for (n2 = 0; n2 < 256; ++n2) {
long double sum = 0.0;
- double lC = aa->is_soap? 0 : lgamma(n1+n2+1) - lgamma(n1+1) - lgamma(n2+1); // \binom{n1+n2}{n1}
+ double lC = aa->errmod == BAM_ERRMOD_SOAP? 0 : lgamma(n1+n2+1) - lgamma(n1+1) - lgamma(n2+1);
for (k = 1; k <= aa->n_hap - 1; ++k) {
double pk = 1.0 / k / sum_harmo;
double log1 = log((double)k/aa->n_hap);
@@ -62,6 +64,7 @@ static void cal_coef(bam_maqcns_t *aa)
long double sum_a[257], b[256], q_c[256], tmp[256], fk2[256];
double *lC;
+ if (aa->errmod == BAM_ERRMOD_MAQ2) return; // no need to do the following
// aa->lhet will be allocated and initialized
free(aa->fk); free(aa->coef);
aa->coef = 0;
@@ -71,7 +74,7 @@ static void cal_coef(bam_maqcns_t *aa)
aa->fk[n] = pow(aa->theta, n) * (1.0 - aa->eta) + aa->eta;
fk2[n] = aa->fk[n>>1]; // this is an approximation, assuming reads equally likely come from both strands
}
- if (aa->is_soap) return;
+ if (aa->errmod == BAM_ERRMOD_SOAP) return;
aa->coef = (double*)calloc(256*256*64, sizeof(double));
lC = (double*)calloc(256 * 256, sizeof(double));
for (n = 1; n != 256; ++n)
@@ -116,19 +119,21 @@ bam_maqcns_t *bam_maqcns_init()
void bam_maqcns_prepare(bam_maqcns_t *bm)
{
+ if (bm->errmod == BAM_ERRMOD_MAQ2) bm->aux->em = errmod_init(1. - bm->theta);
cal_coef(bm); cal_het(bm);
}
void bam_maqcns_destroy(bam_maqcns_t *bm)
{
if (bm == 0) return;
- free(bm->lhet); free(bm->fk); free(bm->coef); free(bm->aux->info);
+ free(bm->lhet); free(bm->fk); free(bm->coef); free(bm->aux->info); free(bm->aux->info16);
+ if (bm->aux->em) errmod_destroy(bm->aux->em);
free(bm->aux); free(bm);
}
glf1_t *bam_maqcns_glfgen(int _n, const bam_pileup1_t *pl, uint8_t ref_base, bam_maqcns_t *bm)
{
- glf_call_aux_t *b;
+ glf_call_aux_t *b = 0;
int i, j, k, w[8], c, n;
glf1_t *g = (glf1_t*)calloc(1, sizeof(glf1_t));
float p[16], min_p = 1e30;
@@ -142,28 +147,38 @@ glf1_t *bam_maqcns_glfgen(int _n, const bam_pileup1_t *pl, uint8_t ref_base, bam
bm->aux->max = _n;
kroundup32(bm->aux->max);
bm->aux->info = (uint32_t*)realloc(bm->aux->info, 4 * bm->aux->max);
+ bm->aux->info16 = (uint16_t*)realloc(bm->aux->info16, 2 * bm->aux->max);
}
- for (i = n = 0; i < _n; ++i) {
+ for (i = n = 0, rms = 0; i < _n; ++i) {
const bam_pileup1_t *p = pl + i;
uint32_t q, x = 0, qq;
+ uint16_t y = 0;
if (p->is_del || (p->b->core.flag&BAM_FUNMAP)) continue;
q = (uint32_t)bam1_qual(p->b)[p->qpos];
x |= (uint32_t)bam1_strand(p->b) << 18 | q << 8 | p->b->core.qual;
+ y |= bam1_strand(p->b)<<4;
if (p->b->core.qual < q) q = p->b->core.qual;
+ c = p->b->core.qual < bm->cap_mapQ? p->b->core.qual : bm->cap_mapQ;
+ rms += c * c;
x |= q << 24;
+ y |= q << 5;
qq = bam1_seqi(bam1_seq(p->b), p->qpos);
q = bam_nt16_nt4_table[qq? qq : ref_base];
- if (!p->is_del && q < 4) x |= 1 << 21 | q << 16;
+ if (!p->is_del && q < 4) x |= 1 << 21 | q << 16, y |= q;
+ bm->aux->info16[n] = y;
bm->aux->info[n++] = x;
}
+ rms = (uint8_t)(sqrt((double)rms / n) + .499);
+ if (bm->errmod == BAM_ERRMOD_MAQ2) {
+ errmod_cal(bm->aux->em, n, 4, bm->aux->info16, p);
+ goto goto_glf;
+ }
ks_introsort(uint32_t, n, bm->aux->info);
// generate esum and fsum
b = (glf_call_aux_t*)calloc(1, sizeof(glf_call_aux_t));
for (k = 0; k != 8; ++k) w[k] = 0;
- rms = 0;
for (j = n - 1; j >= 0; --j) { // calculate esum and fsum
uint32_t info = bm->aux->info[j];
- int tmp;
if (info>>24 < 4 && (info>>8&0x3f) != 0) info = 4<<24 | (info&0xffffff);
k = info>>16&7;
if (info>>24 > 0) {
@@ -172,17 +187,14 @@ glf1_t *bam_maqcns_glfgen(int _n, const bam_pileup1_t *pl, uint8_t ref_base, bam
if (w[k] < 0xff) ++w[k];
++b->c[k&3];
}
- tmp = (int)(info&0xff) < bm->cap_mapQ? (int)(info&0xff) : bm->cap_mapQ;
- rms += tmp * tmp;
}
- b->rms_mapQ = (uint8_t)(sqrt((double)rms / n) + .499);
// rescale ->c[]
for (j = c = 0; j != 4; ++j) c += b->c[j];
if (c > 255) {
for (j = 0; j != 4; ++j) b->c[j] = (int)(254.0 * b->c[j] / c + 0.5);
for (j = c = 0; j != 4; ++j) c += b->c[j];
}
- if (!bm->is_soap) {
+ if (bm->errmod == BAM_ERRMOD_MAQ) {
// generate likelihood
for (j = 0; j != 4; ++j) {
// homozygous
@@ -234,7 +246,7 @@ glf1_t *bam_maqcns_glfgen(int _n, const bam_pileup1_t *pl, uint8_t ref_base, bam
if (max1 > max2 && (min_k != max_k || min1 + 1.0 > min2))
p[max_k<<2|max_k] = min1 > 1.0? min1 - 1.0 : 0.0;
}
- } else { // apply the SOAP model
+ } else if (bm->errmod == BAM_ERRMOD_SOAP) { // apply the SOAP model
// generate likelihood
for (j = 0; j != 4; ++j) {
float tmp;
@@ -251,8 +263,9 @@ glf1_t *bam_maqcns_glfgen(int _n, const bam_pileup1_t *pl, uint8_t ref_base, bam
}
}
+goto_glf:
// convert necessary information to glf1_t
- g->ref_base = ref_base; g->max_mapQ = b->rms_mapQ;
+ g->ref_base = ref_base; g->max_mapQ = rms;
g->depth = n > 16777215? 16777215 : n;
for (j = 0; j != 4; ++j)
for (k = j; k < 4; ++k)
diff --git a/bam_maqcns.h b/bam_maqcns.h
index 6cc5355..28733a5 100644
--- a/bam_maqcns.h
+++ b/bam_maqcns.h
@@ -3,11 +3,15 @@
#include "glf.h"
+#define BAM_ERRMOD_MAQ2 0
+#define BAM_ERRMOD_MAQ 1
+#define BAM_ERRMOD_SOAP 2
+
struct __bmc_aux_t;
typedef struct {
float het_rate, theta;
- int n_hap, cap_mapQ, is_soap;
+ int n_hap, cap_mapQ, errmod;
float eta, q_r;
double *fk, *coef;
diff --git a/bam_plcmd.c b/bam_plcmd.c
index 0cf081d..aef10d3 100644
--- a/bam_plcmd.c
+++ b/bam_plcmd.c
@@ -343,12 +343,12 @@ int bam_pileup(int argc, char *argv[])
d->max_depth = 0;
d->tid = -1; d->mask = BAM_DEF_MASK;
d->c = bam_maqcns_init();
- d->c->is_soap = 1; // change the default model
+ d->c->errmod = BAM_ERRMOD_MAQ2; // change the default model
d->ido = bam_maqindel_opt_init();
while ((c = getopt(argc, argv, "st:f:cT:N:r:l:d:im:gI:G:vM:S2aR:PA")) >= 0) {
switch (c) {
- case 'a': d->c->is_soap = 1; break;
- case 'A': d->c->is_soap = 0; break;
+ case 'a': d->c->errmod = BAM_ERRMOD_SOAP; break;
+ case 'A': d->c->errmod = BAM_ERRMOD_MAQ; break;
case 's': d->format |= BAM_PLF_SIMPLE; break;
case 't': fn_list = strdup(optarg); break;
case 'l': fn_pos = strdup(optarg); break;
diff --git a/errmod.c b/errmod.c
new file mode 100644
index 0000000..fba9a8d
--- /dev/null
+++ b/errmod.c
@@ -0,0 +1,130 @@
+#include <math.h>
+#include "errmod.h"
+#include "ksort.h"
+KSORT_INIT_GENERIC(uint16_t)
+
+typedef struct __errmod_coef_t {
+ double *fk, *beta, *lhet;
+} errmod_coef_t;
+
+typedef struct {
+ double fsum[16], bsum[16];
+ uint32_t c[16];
+} call_aux_t;
+
+static errmod_coef_t *cal_coef(double depcorr, double eta)
+{
+ int k, n, q;
+ long double sum, sum1;
+ double *lC;
+ errmod_coef_t *ec;
+
+ ec = calloc(1, sizeof(errmod_coef_t));
+ // initialize ->fk
+ ec->fk = (double*)calloc(256, sizeof(double));
+ ec->fk[0] = 1.0;
+ for (n = 1; n != 256; ++n)
+ ec->fk[n] = pow(1. - depcorr, n) * (1.0 - eta) + eta;
+ // initialize ->coef
+ ec->beta = (double*)calloc(256 * 256 * 64, sizeof(double));
+ lC = (double*)calloc(256 * 256, sizeof(double));
+ for (n = 1; n != 256; ++n) {
+ double lgn = lgamma(n+1);
+ for (k = 1; k <= n; ++k)
+ lC[n<<8|k] = lgn - lgamma(k+1) - lgamma(n-k+1);
+ }
+ for (q = 1; q != 64; ++q) {
+ double e = pow(10.0, -q/10.0);
+ double le = log(e);
+ double le1 = log(1.0 - e);
+ for (n = 1; n <= 255; ++n) {
+ double *beta = ec->beta + (q<<16|n<<8);
+ sum1 = sum = 0.0;
+ for (k = n; k >= 0; --k, sum1 = sum) {
+ sum = sum1 + expl(lC[n<<8|k] + k*le + (n-k)*le1);
+ beta[k] = -10. / M_LN10 * logl(sum1 / sum);
+ }
+ }
+ }
+ // initialize ->lhet
+ ec->lhet = (double*)calloc(256 * 256, sizeof(double));
+ for (n = 0; n < 256; ++n)
+ for (k = 0; k < 256; ++k)
+ ec->lhet[n<<8|k] = lC[n<<8|k] - M_LN2 * n;
+ free(lC);
+ return ec;
+}
+
+errmod_t *errmod_init(float depcorr)
+{
+ errmod_t *em;
+ em = (errmod_t*)calloc(1, sizeof(errmod_t));
+ em->depcorr = depcorr;
+ em->coef = cal_coef(depcorr, 0.03);
+ return em;
+}
+
+void errmod_destroy(errmod_t *em)
+{
+ if (em == 0) return;
+ free(em->coef->lhet); free(em->coef->fk); free(em->coef->beta);
+ free(em->coef); free(em);
+}
+// qual:6, strand:1, base:4
+int errmod_cal(const errmod_t *em, int n, int m, uint16_t *bases, float *q)
+{
+ call_aux_t aux;
+ int i, j, k, w[32];
+
+ if (m > m) return -1;
+ memset(q, 0, m * m * sizeof(float));
+ if (n == 0) return 0;
+ // calculate aux.esum and aux.fsum
+ if (n > 255) { // then sample 255 bases
+ ks_shuffle(uint16_t, n, bases);
+ n = 255;
+ }
+ ks_introsort(uint16_t, n, bases);
+ memset(w, 0, 32 * sizeof(int));
+ memset(&aux, 0, sizeof(call_aux_t));
+ for (j = n - 1; j >= 0; --j) { // calculate esum and fsum
+ uint16_t b = bases[j];
+ int q = b>>5 < 4? 4 : b>>5;
+ if (q > 63) q = 63;
+ k = b&0x1f;
+ aux.fsum[k&0xf] += em->coef->fk[w[k]];
+ aux.bsum[k&0xf] += em->coef->fk[w[k]] * em->coef->beta[q<<16|n<<8|aux.c[k&0xf]];
+ ++aux.c[k&0xf];
+ ++w[k];
+ }
+ // generate likelihood
+ for (j = 0; j != m; ++j) {
+ float tmp1, tmp3;
+ int tmp2, bar_e;
+ // homozygous
+ for (k = 0, tmp1 = tmp3 = 0.0, tmp2 = 0; k != m; ++k) {
+ if (k == j) continue;
+ tmp1 += aux.bsum[k]; tmp2 += aux.c[k]; tmp3 += aux.fsum[k];
+ }
+ if (tmp2) {
+ bar_e = (int)(tmp1 / tmp3 + 0.499);
+ if (bar_e > 63) bar_e = 63;
+ q[j*m+j] = tmp1;
+ }
+ // heterozygous
+ for (k = j + 1; k < m; ++k) {
+ int cjk = aux.c[j] + aux.c[k];
+ for (i = 0, tmp2 = 0, tmp1 = tmp3 = 0.0; i < m; ++i) {
+ if (i == j || i == k) continue;
+ tmp1 += aux.bsum[i]; tmp2 += aux.c[i]; tmp3 += aux.fsum[i];
+ }
+ if (tmp2) {
+ bar_e = (int)(tmp1 / tmp3 + 0.499);
+ if (bar_e > 63) bar_e = 63;
+ q[j*m+k] = q[k*m+j] = -4.343 * em->coef->lhet[cjk<<8|aux.c[k]] + tmp1;
+ } else q[j*m+k] = q[k*m+j] = -4.343 * em->coef->lhet[cjk<<8|aux.c[k]]; // all the bases are either j or k
+ }
+ for (k = 0; k != m; ++k) if (q[j*m+k] < 0.0) q[j*m+k] = 0.0;
+ }
+ return 0;
+}
diff --git a/errmod.h b/errmod.h
new file mode 100644
index 0000000..e3e9a90
--- /dev/null
+++ b/errmod.h
@@ -0,0 +1,17 @@
+#ifndef ERRMOD_H
+#define ERRMOD_H
+
+#include <stdint.h>
+
+struct __errmod_coef_t;
+
+typedef struct {
+ double depcorr;
+ struct __errmod_coef_t *coef;
+} errmod_t;
+
+errmod_t *errmod_init(float depcorr);
+void errmod_destroy(errmod_t *em);
+int errmod_cal(const errmod_t *em, int n, int m, uint16_t *bases, float *q);
+
+#endif
diff --git a/ksort.h b/ksort.h
index 16a03fd..fa850ab 100644
--- a/ksort.h
+++ b/ksort.h
@@ -250,6 +250,15 @@ typedef struct {
if (hh <= k) low = ll; \
if (hh >= k) high = hh - 1; \
} \
+ } \
+ void ks_shuffle_##name(size_t n, type_t a[]) \
+ { \
+ int i, j; \
+ for (i = n; i > 1; --i) { \
+ type_t tmp; \
+ j = (int)(drand48() * i); \
+ tmp = a[j]; a[j] = a[i-1]; a[i-1] = tmp; \
+ } \
}
#define ks_mergesort(name, n, a, t) ks_mergesort_##name(n, a, t)
@@ -259,6 +268,7 @@ typedef struct {
#define ks_heapmake(name, n, a) ks_heapmake_##name(n, a)
#define ks_heapadjust(name, i, n, a) ks_heapadjust_##name(i, n, a)
#define ks_ksmall(name, n, a, k) ks_ksmall_##name(n, a, k)
+#define ks_shuffle(name, n, a) ks_shuffle_##name(n, a)
#define ks_lt_generic(a, b) ((a) < (b))
#define ks_lt_str(a, b) (strcmp((a), (b)) < 0)
--
2.39.5