From d14782f26230244039bf38f6384a7a9106b0daed Mon Sep 17 00:00:00 2001 From: Heng Li Date: Mon, 23 Aug 2010 18:06:07 +0000 Subject: [PATCH] * samtools-0.1.8-12 (r692) * group data by samples in "mpileup -g" --- Makefile | 2 +- bam_plcmd.c | 73 +++++++++++++++++++++++++++++-------- bamtk.c | 2 +- examples/Makefile | 2 +- sample.c | 92 +++++++++++++++++++++++++++++++++++++++++++++++ sample.h | 17 +++++++++ 6 files changed, 170 insertions(+), 18 deletions(-) create mode 100644 sample.c create mode 100644 sample.h diff --git a/Makefile b/Makefile index 6a15c28..3ff3272 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 bam2bcf.o errmod.o + bamtk.o kaln.o bam2bcf.o errmod.o sample.o PROG= samtools INCLUDES= -I. SUBDIRS= . bcftools misc diff --git a/bam_plcmd.c b/bam_plcmd.c index 80799e3..a23477f 100644 --- a/bam_plcmd.c +++ b/bam_plcmd.c @@ -5,7 +5,6 @@ #include "sam.h" #include "faidx.h" #include "bam_maqcns.h" -#include "bam2bcf.h" #include "khash.h" #include "glf.h" #include "kstring.h" @@ -449,6 +448,10 @@ int bam_pileup(int argc, char *argv[]) * mpileup * ***********/ +#include +#include "bam2bcf.h" +#include "sample.h" + #define MPLP_GLF 0x10 #define MPLP_NO_COMP 0x20 @@ -466,6 +469,12 @@ typedef struct { int min_mq; } mplp_aux_t; +typedef struct { + int n; + int *n_plp, *m_plp; + bam_pileup1_t **plp; +} mplp_pileup_t; + static int mplp_func(void *data, bam1_t *b) { mplp_aux_t *ma = (mplp_aux_t*)data; @@ -476,6 +485,29 @@ static int mplp_func(void *data, bam1_t *b) return ret; } +static void group_smpl(mplp_pileup_t *m, bam_sample_t *sm, kstring_t *buf, + int n, char *const*fn, int *n_plp, const bam_pileup1_t **plp) +{ + int i, j; + memset(m->n_plp, 0, m->n * sizeof(int)); + for (i = 0; i < n; ++i) { + for (j = 0; j < n_plp[i]; ++j) { + const bam_pileup1_t *p = plp[i] + j; + uint8_t *q; + int id = -1; + q = bam_aux_get(p->b, "RG"); + if (q) id = bam_smpl_rg2smid(sm, fn[i], (char*)q+1, buf); + if (id < 0) id = bam_smpl_rg2smid(sm, fn[i], 0, buf); + assert(id >= 0 && id < m->n); + if (m->n_plp[id] == m->m_plp[id]) { + m->m_plp[id] = m->m_plp[id]? m->m_plp[id]<<1 : 8; + m->plp[id] = realloc(m->plp[id], sizeof(bam_pileup1_t) * m->m_plp[id]); + } + m->plp[id][m->n_plp[id]++] = *p; + } + } +} + static int mpileup(mplp_conf_t *conf, int n, char **fn) { mplp_aux_t **data; @@ -492,10 +524,17 @@ static int mpileup(mplp_conf_t *conf, int n, char **fn) bcf_t *bp = 0; bcf_hdr_t *bh = 0; + bam_sample_t *sm = 0; + kstring_t buf; + mplp_pileup_t gplp; + + memset(&gplp, 0, sizeof(mplp_pileup_t)); + memset(&buf, 0, sizeof(kstring_t)); memset(&bc, 0, sizeof(bcf_call_t)); data = calloc(n, sizeof(void*)); plp = calloc(n, sizeof(void*)); n_plp = calloc(n, sizeof(int*)); + sm = bam_smpl_init(); // read the header and initialize data for (i = 0; i < n; ++i) { @@ -504,6 +543,7 @@ static int mpileup(mplp_conf_t *conf, int n, char **fn) data[i]->min_mq = conf->min_mq; data[i]->fp = bam_open(fn[i], "r"); h_tmp = bam_header_read(data[i]->fp); + bam_smpl_add(sm, fn[i], h_tmp->text); if (conf->reg) { int beg, end; bam_index_t *idx; @@ -526,6 +566,12 @@ static int mpileup(mplp_conf_t *conf, int n, char **fn) bam_header_destroy(h_tmp); } } + gplp.n = sm->n; + gplp.n_plp = calloc(sm->n, sizeof(int)); + gplp.m_plp = calloc(sm->n, sizeof(int)); + gplp.plp = calloc(sm->n, sizeof(void*)); + + fprintf(stderr, "[%s] %d samples in %d input files\n", __func__, sm->n, n); if (conf->fn_pos) hash = load_pos(conf->fn_pos, h); // write the VCF header if (conf->flag & MPLP_GLF) { @@ -541,12 +587,8 @@ static int mpileup(mplp_conf_t *conf, int n, char **fn) bh->name = malloc(s.l); memcpy(bh->name, s.s, s.l); s.l = 0; - for (i = 0; i < n; ++i) { - const char *p; - if ((p = strstr(fn[i], ".bam")) != 0) - kputsn(fn[i], p - fn[i], &s); - else kputs(fn[i], &s); - kputc('\0', &s); + for (i = 0; i < sm->n; ++i) { + kputs(sm->smpl[i], &s); kputc('\0', &s); } bh->l_smpl = s.l; bh->sname = malloc(s.l); @@ -555,11 +597,8 @@ static int mpileup(mplp_conf_t *conf, int n, char **fn) free(s.s); bcf_hdr_sync(bh); bcf_hdr_write(bp, bh); - } - // mpileup - if (conf->flag & MPLP_GLF) { bca = bcf_call_init(-1., conf->min_baseQ); - bcr = calloc(n, sizeof(bcf_callret1_t)); + bcr = calloc(sm->n, sizeof(bcf_callret1_t)); } ref_tid = -1; ref = 0; iter = bam_mplp_init(n, mplp_func, (void**)data); @@ -578,14 +617,14 @@ static int mpileup(mplp_conf_t *conf, int n, char **fn) if (conf->flag & MPLP_GLF) { int _ref0, ref16; bcf1_t *b = calloc(1, sizeof(bcf1_t)); + group_smpl(&gplp, sm, &buf, n, fn, n_plp, plp); _ref0 = (ref && pos < ref_len)? ref[pos] : 'N'; ref16 = bam_nt16_table[_ref0]; - for (i = 0; i < n; ++i) - bcf_call_glfgen(n_plp[i], plp[i], ref16, bca, bcr + i); - bcf_call_combine(n, bcr, ref16, &bc); + for (i = 0; i < gplp.n; ++i) + bcf_call_glfgen(gplp.n_plp[i], gplp.plp[i], ref16, bca, bcr + i); + bcf_call_combine(gplp.n, bcr, ref16, &bc); bcf_call2bcf(tid, pos, &bc, b); bcf_write(bp, bh, b); - //fprintf(stderr, "%d,%d,%d\n", b->tid, b->pos, b->l_str); bcf_destroy(b); } else { printf("%s\t%d\t%c", h->target_name[tid], pos + 1, (ref && pos < ref_len)? ref[pos] : 'N'); @@ -608,7 +647,11 @@ static int mpileup(mplp_conf_t *conf, int n, char **fn) putchar('\n'); } } + bcf_close(bp); + bam_smpl_destroy(sm); free(buf.s); + for (i = 0; i < gplp.n; ++i) free(gplp.plp[i]); + free(gplp.plp); free(gplp.n_plp); free(gplp.m_plp); if (hash) { // free the hash table khint_t k; for (k = kh_begin(hash); k < kh_end(hash); ++k) diff --git a/bamtk.c b/bamtk.c index 480f746..d4b40ea 100644 --- a/bamtk.c +++ b/bamtk.c @@ -9,7 +9,7 @@ #endif #ifndef PACKAGE_VERSION -#define PACKAGE_VERSION "0.1.8-11 (r672)" +#define PACKAGE_VERSION "0.1.8-12 (r693)" #endif int bam_taf2baf(int argc, char *argv[]); diff --git a/examples/Makefile b/examples/Makefile index 96a4e7e..2733a57 100644 --- a/examples/Makefile +++ b/examples/Makefile @@ -19,7 +19,7 @@ ex1a.bam:ex1.bam ex1b.bam:ex1.bam ../samtools view -h ex1.bam | awk 'BEGIN{FS=OFS="\t"}{if(/^@/)print;else{$$1=$$1"b";print}}' | ../samtools view -bS - > $@ ex1f.rg: - (echo "@RG ID:ex1 LB:ex1"; echo "@RG ID:ex1a LB:ex1"; echo "@RG ID:ex1b LB:ex1b") > $@ + (echo "@RG ID:ex1 LB:ex1 SM:ex1"; echo "@RG ID:ex1a LB:ex1 SM:ex1"; echo "@RG ID:ex1b LB:ex1b SM:ex1b") > $@ ex1f.bam:ex1.bam ex1a.bam ex1b.bam ex1f.rg ../samtools merge -rh ex1f.rg $@ ex1.bam ex1a.bam ex1b.bam ex1f-rmduppe.bam:ex1f.bam diff --git a/sample.c b/sample.c new file mode 100644 index 0000000..95bec68 --- /dev/null +++ b/sample.c @@ -0,0 +1,92 @@ +#include +#include +#include "sample.h" +#include "khash.h" +KHASH_MAP_INIT_STR(sm, int) + +bam_sample_t *bam_smpl_init(void) +{ + bam_sample_t *s; + s = calloc(1, sizeof(bam_sample_t)); + s->rg2smid = kh_init(sm); + return s; +} + +void bam_smpl_destroy(bam_sample_t *sm) +{ + int i; + khint_t k; + khash_t(sm) *rg2smid = (khash_t(sm)*)sm->rg2smid; + if (sm == 0) return; + for (i = 0; i < sm->n; ++i) free(sm->smpl[i]); + free(sm->smpl); + for (k = kh_begin(rg2smid); k != kh_end(rg2smid); ++k) + if (kh_exist(rg2smid, k)) free((char*)kh_key(rg2smid, k)); + kh_destroy(sm, sm->rg2smid); + free(sm); +} + +static void add_pair(bam_sample_t *sm, khash_t(sm) *sm2id, const char *key, const char *val) +{ + khint_t k_rg, k_sm; + int ret; + khash_t(sm) *rg2smid = (khash_t(sm)*)sm->rg2smid; + k_rg = kh_get(sm, rg2smid, key); + if (k_rg != kh_end(rg2smid)) return; // duplicated @RG-ID + k_rg = kh_put(sm, rg2smid, strdup(key), &ret); + k_sm = kh_get(sm, sm2id, val); + if (k_sm == kh_end(sm2id)) { // absent + if (sm->n == sm->m) { + sm->m = sm->m? sm->m<<1 : 1; + sm->smpl = realloc(sm->smpl, sizeof(void*) * sm->m); + } + sm->smpl[sm->n] = strdup(val); + k_sm = kh_put(sm, sm2id, sm->smpl[sm->n], &ret); + kh_val(sm2id, k_sm) = sm->n++; + } + kh_val(rg2smid, k_rg) = kh_val(sm2id, k_sm); +} + +int bam_smpl_add(bam_sample_t *sm, const char *fn, const char *txt) +{ + const char *p = txt, *q, *r; + kstring_t buf; + int n = 0; + khash_t(sm) *sm2id; + sm2id = kh_init(sm); + memset(&buf, 0, sizeof(kstring_t)); + while ((q = strstr(p, "@RG")) != 0) { + p = q + 3; + r = q = 0; + if ((q = strstr(p, "\tID:")) != 0) q += 4; + if ((r = strstr(p, "\tSM:")) != 0) r += 4; + if (r && q) { + char *u, *v; + int oq, or; + for (u = (char*)q; *u && *u != '\t' && *u != '\n'; ++u); + for (v = (char*)r; *v && *v != '\t' && *v != '\n'; ++v); + oq = *u; or = *v; *u = *v = '\0'; + buf.l = 0; kputs(fn, &buf); kputc('/', &buf); kputs(q, &buf); + add_pair(sm, sm2id, buf.s, r); + *u = oq; *v = or; + } else break; + p = q > r? q : r; + ++n; + } + if (n == 0) add_pair(sm, sm2id, fn, fn); + free(buf.s); + kh_destroy(sm, sm2id); + return 0; +} + +int bam_smpl_rg2smid(const bam_sample_t *sm, const char *fn, const char *rg, kstring_t *str) +{ + khint_t k; + khash_t(sm) *rg2smid = (khash_t(sm)*)sm->rg2smid; + if (rg) { + str->l = 0; + kputs(fn, str); kputc('/', str); kputs(rg, str); + k = kh_get(sm, rg2smid, str->s); + } else k = kh_get(sm, rg2smid, fn); + return k == kh_end(rg2smid)? -1 : kh_val(rg2smid, k); +} diff --git a/sample.h b/sample.h new file mode 100644 index 0000000..149f380 --- /dev/null +++ b/sample.h @@ -0,0 +1,17 @@ +#ifndef BAM_SAMPLE_H +#define BAM_SAMPLE_H + +#include "kstring.h" + +typedef struct { + int n, m; + char **smpl; + void *rg2smid; +} bam_sample_t; + +bam_sample_t *bam_smpl_init(void); +int bam_smpl_add(bam_sample_t *sm, const char *abs, const char *txt); +int bam_smpl_rg2smid(const bam_sample_t *sm, const char *fn, const char *rg, kstring_t *str); +void bam_smpl_destroy(bam_sample_t *sm); + +#endif -- 2.39.2