From 921db5eb2dfecd81e71c42934d0482efb869ab80 Mon Sep 17 00:00:00 2001 From: Heng Li Date: Wed, 14 Jul 2010 17:33:51 +0000 Subject: [PATCH] * added mutli-sample framework. It is not working, yet. * improved the mpileup interface --- Makefile | 2 +- bam_mcns.c | 127 ++++++++++++++++++++++++++++++++++++++++++++++++++++ bam_mcns.h | 22 +++++++++ bam_plcmd.c | 91 ++++++++++++++++++++++++++++++------- 4 files changed, 224 insertions(+), 18 deletions(-) create mode 100644 bam_mcns.c create mode 100644 bam_mcns.h diff --git a/Makefile b/Makefile index 35d578f..88d8217 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 + bamtk.o kaln.o bam_mcns.o PROG= samtools INCLUDES= SUBDIRS= . misc diff --git a/bam_mcns.c b/bam_mcns.c new file mode 100644 index 0000000..845dfcd --- /dev/null +++ b/bam_mcns.c @@ -0,0 +1,127 @@ +#include +#include +#include "bam_mcns.h" + +#define MC_MIN_QUAL 20 +#define MC_MAX_SUMQ 3000 +#define MC_MAX_SUMQP 1e-300 + +struct __mc_aux_t { + int n, N; + int ref, alt; + double *q2p, *pdg; // pdg -> P(D|g) + int *qsum, *bcnt; +}; + +mc_aux_t *mc_init(int n) // FIXME: assuming diploid +{ + mc_aux_t *ma; + int i; + ma = calloc(1, sizeof(mc_aux_t)); + ma->n = n; ma->N = 2 * n; + ma->q2p = calloc(MC_MAX_SUMQ + 1, sizeof(double)); + ma->qsum = calloc(4 * ma->n, sizeof(int)); + ma->bcnt = calloc(4 * ma->n, sizeof(int)); + ma->pdg = calloc(3 * ma->n, sizeof(double)); + for (i = 0; i <= MC_MAX_SUMQ; ++i) + ma->q2p[i] = pow(10., -i / 10.); + return ma; +} + +void mc_destroy(mc_aux_t *ma) +{ + if (ma) { + free(ma->qsum); free(ma->bcnt); + free(ma->q2p); free(ma->pdg); + free(ma); + } +} + +static void sum_err(int *n, const bam_pileup1_t **plp, mc_aux_t *ma) +{ + int i, j; + memset(ma->qsum, 0, sizeof(int) * 4 * ma->n); + memset(ma->bcnt, 0, sizeof(int) * 4 * ma->n); + for (j = 0; j < ma->n; ++j) { + int *qsum = ma->qsum + j * 4; + int *bcnt = ma->bcnt + j * 4; + for (i = 0; i < n[j]; ++i) { + const bam_pileup1_t *p = plp[j] + i; + int q, b; + if (p->is_del || (p->b->core.flag&BAM_FUNMAP)) continue; + q = bam1_qual(p->b)[p->qpos]; + if (p->b->core.qual < q) q = p->b->core.qual; + if (q < MC_MIN_QUAL) continue; // small qual + b = bam_nt16_nt4_table[(int)bam1_seqi(bam1_seq(p->b), p->qpos)]; + if (b > 3) continue; // N + qsum[b] += q; + ++bcnt[b]; + } + } +} + +static void set_allele(int ref, mc_aux_t *ma) +{ + int i, j, sum[4], tmp; + sum[0] = sum[1] = sum[2] = sum[3] = 0; + for (i = 0; i < ma->n; ++i) + for (j = 0; j < 4; ++j) + sum[j] += ma->qsum[i * 4 + j]; + for (j = 0; j < 4; ++j) sum[j] = sum[j]<<2 | j; + for (i = 1; i < 4; ++i) // insertion sort + for (j = i; j > 0 && sum[j] < sum[j-1]; --j) + tmp = sum[j], sum[j] = sum[j-1], sum[j-1] = tmp; + ma->ref = sum[3]&3; ma->alt = sum[2]&3; + if (ref == ma->alt) tmp = ma->ref, ma->ref = ma->alt, ma->alt = tmp; + // note that ma->ref might not be ref in case of triallele +} + +static void cal_pdg(mc_aux_t *ma) +{ + int i, j; + for (j = 0; j < ma->n; ++j) { + int pi[3], *qsum, *bcnt; + double *pdg = ma->pdg + j * 3; + qsum = ma->qsum + j * 4; + bcnt = ma->bcnt + j * 4; + pi[1] = 3 * (bcnt[ma->ref] + bcnt[ma->alt]); + pi[0] = qsum[ma->alt]; + pi[2] = qsum[ma->ref]; + for (i = 0; i < 3; ++i) + pdg[i] = pi[i] < MC_MAX_SUMQ? MC_MAX_SUMQP : ma->q2p[pi[i]]; + } +} + +double mc_freq0(int ref, int *n, const bam_pileup1_t **plp, mc_aux_t *ma, int *_ref, int *alt) +{ + int i, acnt[4], j; + double f0; + sum_err(n, plp, ma); + set_allele(ref, ma); + cal_pdg(ma); + acnt[0] = acnt[1] = acnt[2] = acnt[3] = 0; + for (i = 0; i < ma->n; ++i) + for (j = 0; j < 4; ++j) + acnt[j] += ma->bcnt[i * 4 + j]; + f0 = acnt[ma->ref] + acnt[ma->alt] == 0? -1. + : (double)acnt[ref] / (acnt[ma->ref] + acnt[ma->alt]); + *_ref = ma->ref; *alt = ma->alt; + return f0; +} + +double mc_freq_iter(double f0, mc_aux_t *ma) +{ + double f, f3[3]; + int i, j; + f3[0] = f0*f0; f3[1] = 2.*f0*(1.-f0); f3[2] = (1.-f0)*(1.-f0); + for (i = 0, f = 0.; i < ma->n; ++i) { + double up, dn, *pdg; + pdg = ma->pdg + i * 3; + for (j = 1, up = 0.; j < 3; ++j) + up += j * pdg[j] * f3[j]; + for (j = 0, dn = 0.; j < 3; ++j) + dn += pdg[j] * f3[j]; + f += up / dn; + } + return f; +} diff --git a/bam_mcns.h b/bam_mcns.h new file mode 100644 index 0000000..ac1e41c --- /dev/null +++ b/bam_mcns.h @@ -0,0 +1,22 @@ +#ifndef BAM_MCNS_H +#define BAM_MCNS_H + +#include "bam.h" + +struct __mc_aux_t; +typedef struct __mc_aux_t mc_aux_t; + +#ifdef __cplusplus +extern "C" { +#endif + + mc_aux_t *mc_init(int n); + void mc_destroy(mc_aux_t *ma); + double mc_freq0(int ref, int *n, const bam_pileup1_t **plp, mc_aux_t *ma, int *_ref, int *alt); + double mc_freq_iter(double f0, mc_aux_t *ma); + +#ifdef __cplusplus +} +#endif + +#endif diff --git a/bam_plcmd.c b/bam_plcmd.c index c665a92..e4de833 100644 --- a/bam_plcmd.c +++ b/bam_plcmd.c @@ -5,6 +5,7 @@ #include "sam.h" #include "faidx.h" #include "bam_maqcns.h" +#include "bam_mcns.h" #include "khash.h" #include "glf.h" #include "kstring.h" @@ -449,6 +450,7 @@ int bam_pileup(int argc, char *argv[]) ***********/ typedef struct { + int vcf, max_mq, min_mq; char *reg, *fn_pos; faidx_t *fai; kh_64_t *hash; @@ -457,18 +459,23 @@ typedef struct { typedef struct { bamFile fp; bam_iter_t iter; + int min_mq; } mplp_aux_t; static int mplp_func(void *data, bam1_t *b) { mplp_aux_t *ma = (mplp_aux_t*)data; - if (ma->iter) return bam_iter_read(ma->fp, ma->iter, b); - return bam_read1(ma->fp, b); + int ret; + do { + ret = ma->iter? bam_iter_read(ma->fp, ma->iter, b) : bam_read1(ma->fp, b); + } while (b->core.qual < ma->min_mq && ret >= 0); + return ret; } static int mpileup(mplp_conf_t *conf, int n, char **fn) { mplp_aux_t **data; + mc_aux_t *ma = 0; int i, tid, pos, *n_plp, beg0 = 0, end0 = 1u<<29, ref_len, ref_tid; const bam_pileup1_t **plp; bam_mplp_t iter; @@ -483,6 +490,7 @@ static int mpileup(mplp_conf_t *conf, int n, char **fn) for (i = 0; i < n; ++i) { bam_header_t *h_tmp; data[i] = calloc(1, sizeof(mplp_aux_t)); + data[i]->min_mq = conf->min_mq; data[i]->fp = bam_open(fn[i], "r"); h_tmp = bam_header_read(data[i]->fp); if (conf->reg) { @@ -508,7 +516,23 @@ static int mpileup(mplp_conf_t *conf, int n, char **fn) } } if (conf->fn_pos) hash = load_pos(conf->fn_pos, h); + // write the VCF header + if (conf->vcf) { + kstring_t s; + s.l = s.m = 0; s.s = 0; + kputs("#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\tFORMAT", &s); + for (i = 0; i < n; ++i) { + const char *p; + kputc('\t', &s); + if ((p = strstr(fn[i], ".bam")) != 0) + kputsn(fn[i], p - fn[i], &s); + else kputs(fn[i], &s); + } + puts(s.s); + free(s.s); + } // mpileup + if (conf->vcf) ma = mc_init(n); ref_tid = -1; ref = 0; iter = bam_mplp_init(n, mplp_func, (void**)data); while (bam_mplp_auto(iter, &tid, &pos, n_plp, plp) > 0) { @@ -523,24 +547,52 @@ static int mpileup(mplp_conf_t *conf, int n, char **fn) if (conf->fai) ref = fai_fetch(conf->fai, h->target_name[tid], &ref_len); ref_tid = tid; } - 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) { - int j; - printf("\t%d\t", n_plp[i]); - if (n_plp[i] == 0) printf("*\t*"); - else { - for (j = 0; j < n_plp[i]; ++j) - pileup_seq(plp[i] + j, pos, ref_len, ref); - putchar('\t'); + if (conf->vcf) { + double f; + int j, _ref, _alt, _ref0, depth, rms_q; + uint64_t sqr_sum; + _ref0 = (ref && pos < ref_len)? ref[pos] : 'N'; + _ref0 = bam_nt16_nt4_table[bam_nt16_table[_ref0]]; + f = mc_freq0(_ref0, n_plp, plp, ma, &_ref, &_alt); + printf("%s\t%d\t.\t%c\t", h->target_name[tid], pos + 1, "ACGTN"[_ref0]); + if (_ref0 == _ref) putchar("ACGTN"[_alt]); + else printf("%c,%c", "ACGTN"[_ref], "ACGTN"[_alt]); + printf("\t0\tPASS\t"); // FIXME: currently these not available + for (i = depth = 0, sqr_sum = 0; i < n; ++i) { + depth += n_plp[i]; for (j = 0; j < n_plp[i]; ++j) { - const bam_pileup1_t *p = plp[i] + j; - int c = bam1_qual(p->b)[p->qpos] + 33; - if (c > 126) c = 126; - putchar(c); + int q = plp[i][j].b->core.qual; + if (q > conf->max_mq) q = conf->max_mq; + sqr_sum += q * q; + } + } + rms_q = (int)(sqrt((double)sqr_sum / depth) + .499); + printf("DP=%d;MQ=%d;AF=%.3lg", depth, rms_q, 1.-f); // FIXME: not working for triallelic alleles + printf("\tDP"); + // output genotype information; FIXME: to be implmented... + for (i = 0; i < n; ++i) + printf("\t%d", n_plp[i]); + putchar('\n'); + } 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) { + int j; + printf("\t%d\t", n_plp[i]); + if (n_plp[i] == 0) printf("*\t*"); + else { + for (j = 0; j < n_plp[i]; ++j) + pileup_seq(plp[i] + j, pos, ref_len, ref); + putchar('\t'); + for (j = 0; j < n_plp[i]; ++j) { + const bam_pileup1_t *p = plp[i] + j; + int c = bam1_qual(p->b)[p->qpos] + 33; + if (c > 126) c = 126; + putchar(c); + } } } + putchar('\n'); } - putchar('\n'); } if (hash) { // free the hash table khint_t k; @@ -548,6 +600,7 @@ static int mpileup(mplp_conf_t *conf, int n, char **fn) if (kh_exist(hash, k)) free(kh_val(hash, k)); kh_destroy(64, hash); } + mc_destroy(ma); bam_mplp_destroy(iter); bam_header_destroy(h); for (i = 0; i < n; ++i) { @@ -564,7 +617,8 @@ int bam_mpileup(int argc, char *argv[]) int c; mplp_conf_t mplp; memset(&mplp, 0, sizeof(mplp_conf_t)); - while ((c = getopt(argc, argv, "f:r:l:")) >= 0) { + mplp.max_mq = 60; + while ((c = getopt(argc, argv, "f:r:l:VM:q:")) >= 0) { switch (c) { case 'f': mplp.fai = fai_load(optarg); @@ -572,6 +626,9 @@ int bam_mpileup(int argc, char *argv[]) break; case 'r': mplp.reg = strdup(optarg); break; case 'l': mplp.fn_pos = strdup(optarg); break; + case 'V': mplp.vcf = 1; break; + case 'M': mplp.max_mq = atoi(optarg); break; + case 'q': mplp.min_mq = atoi(optarg); break; } } if (argc == 1) { -- 2.39.2