From 4967051aa734ece6ce44f0fb891a936e849c2fe2 Mon Sep 17 00:00:00 2001 From: Heng Li Date: Sat, 12 Jun 2010 22:02:38 +0000 Subject: [PATCH] * samtools-0.1.7-15 (r592) * fixed a few minor memory leaks in the new pileup code * improved the functionality of mpileup --- ChangeLog | 202 +++++++++++++++++++++++++++++++++++++++++++++++++++ bam_aux.c | 2 +- bam_index.c | 5 +- bam_pileup.c | 7 +- bam_plcmd.c | 95 ++++++++++++++++++++---- bamtk.c | 2 +- 6 files changed, 291 insertions(+), 22 deletions(-) diff --git a/ChangeLog b/ChangeLog index 64df707..2e6f535 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,3 +1,205 @@ +------------------------------------------------------------------------ +r591 | lh3lh3 | 2010-06-12 14:09:22 -0400 (Sat, 12 Jun 2010) | 3 lines +Changed paths: + M /trunk/samtools/bam.h + M /trunk/samtools/bam_pileup.c + M /trunk/samtools/bam_plcmd.c + M /trunk/samtools/bamtk.c + + * samtools-0.1.7-14 (r591) + * elementary multi-way pileup. More testing and more functionality to be done. + +------------------------------------------------------------------------ +r590 | lh3lh3 | 2010-06-12 01:00:24 -0400 (Sat, 12 Jun 2010) | 3 lines +Changed paths: + M /trunk/samtools/bam.h + M /trunk/samtools/bam_pileup.c + M /trunk/samtools/bamtk.c + + * samtools-0.1.7-13 (r590) + * added mpileup APIs. No compiling errors, but not tested at all. It is late. + +------------------------------------------------------------------------ +r589 | lh3lh3 | 2010-06-11 22:37:09 -0400 (Fri, 11 Jun 2010) | 3 lines +Changed paths: + M /trunk/samtools/bam.h + M /trunk/samtools/bam_pileup.c + M /trunk/samtools/bamtk.c + + * samtools-0.1.7-12 (r589) + * added iterator-like APIs for pileup + +------------------------------------------------------------------------ +r588 | lh3lh3 | 2010-06-11 17:41:13 -0400 (Fri, 11 Jun 2010) | 3 lines +Changed paths: + M /trunk/samtools/bam_index.c + M /trunk/samtools/bamtk.c + + * samtools-0.1.7-11 (r588) + * ported a few improvements from tabix back to samtools + +------------------------------------------------------------------------ +r587 | lh3lh3 | 2010-06-11 17:33:16 -0400 (Fri, 11 Jun 2010) | 3 lines +Changed paths: + M /trunk/samtools/bam.h + M /trunk/samtools/bam_index.c + M /trunk/samtools/bamtk.c + + * samtools-0.1.7-10 (r587) + * added iterator interface for bam_fetch (ported back from tabix) + +------------------------------------------------------------------------ +r586 | lh3lh3 | 2010-06-11 13:23:53 -0400 (Fri, 11 Jun 2010) | 3 lines +Changed paths: + M /trunk/samtools/Makefile + A /trunk/samtools/bam_reheader.c + M /trunk/samtools/bamtk.c + M /trunk/samtools/bgzf.c + + * samtools-0.1.7-9 (r586) + * added "reheader" to replace the BAM header + +------------------------------------------------------------------------ +r585 | lh3lh3 | 2010-06-11 12:22:06 -0400 (Fri, 11 Jun 2010) | 3 lines +Changed paths: + M /trunk/samtools/bam.c + M /trunk/samtools/bamtk.c + M /trunk/samtools/kstring.h + + * samtools-0.1.7-8 (r585) + * speed up "view" + +------------------------------------------------------------------------ +r584 | lh3lh3 | 2010-06-11 12:00:41 -0400 (Fri, 11 Jun 2010) | 4 lines +Changed paths: + M /trunk/samtools/bam.c + M /trunk/samtools/bamtk.c + M /trunk/samtools/bgzf.c + M /trunk/samtools/bgzf.h + M /trunk/samtools/kstring.h + M /trunk/samtools/misc/wgsim_eval.pl + + * samtools-0.1.7-7 (r584) + * ported tabix BGZF to samtools + * flush BGZF after writing the BAM header and between alignment boundaries + +------------------------------------------------------------------------ +r583 | petulda | 2010-06-11 11:58:20 -0400 (Fri, 11 Jun 2010) | 1 line +Changed paths: + A /trunk/samtools/misc/varfilter.py + +Initial release on behalf of Aylwyn Scally +------------------------------------------------------------------------ +r561 | petulda | 2010-05-07 08:41:56 -0400 (Fri, 07 May 2010) | 1 line +Changed paths: + M /trunk/samtools/samtools.1 + +Added a note about the indels coordinates +------------------------------------------------------------------------ +r551 | petulda | 2010-04-23 09:42:13 -0400 (Fri, 23 Apr 2010) | 1 line +Changed paths: + M /trunk/samtools/misc/sam2vcf.pl + +Added the possibility to print or not to print the reference allele +------------------------------------------------------------------------ +r546 | petulda | 2010-04-15 04:33:55 -0400 (Thu, 15 Apr 2010) | 1 line +Changed paths: + M /trunk/samtools/sam_header.c + +More descriptive message for space separated tags +------------------------------------------------------------------------ +r545 | petulda | 2010-04-14 11:44:50 -0400 (Wed, 14 Apr 2010) | 1 line +Changed paths: + M /trunk/samtools/misc/sam2vcf.pl + +Speedup with -i, no need to query the reference all the time +------------------------------------------------------------------------ +r541 | petulda | 2010-03-15 10:03:51 -0400 (Mon, 15 Mar 2010) | 1 line +Changed paths: + M /trunk/samtools/sam_header.c + +Fixed the order of sequences in the header +------------------------------------------------------------------------ +r540 | petulda | 2010-03-04 06:28:35 -0500 (Thu, 04 Mar 2010) | 1 line +Changed paths: + M /trunk/samtools/misc/sam2vcf.pl + +Added possibility to select indels only and fixed a bug in reporting homozygous indels. +------------------------------------------------------------------------ +r539 | jmarshall | 2010-02-27 06:48:17 -0500 (Sat, 27 Feb 2010) | 4 lines +Changed paths: + M /trunk/samtools/bam.c + +Improve the invalid 'BAM\1' magic number error message, and also print it +when no bytes can be read from the alleged BAM file, e.g., in the common +user error case when a SAM file has accidentally been supplied. + +------------------------------------------------------------------------ +r538 | petulda | 2010-02-26 10:51:40 -0500 (Fri, 26 Feb 2010) | 1 line +Changed paths: + M /trunk/samtools/AUTHORS + M /trunk/samtools/bam.h + M /trunk/samtools/bam_import.c + M /trunk/samtools/sam_header.c + +Improved efficiency of header parsing +------------------------------------------------------------------------ +r537 | lh3lh3 | 2010-02-23 21:08:48 -0500 (Tue, 23 Feb 2010) | 3 lines +Changed paths: + M /trunk/samtools/misc/export2sam.pl + +Updated export2sam.pl by Chris Saunders from Illumina. + + +------------------------------------------------------------------------ +r536 | petulda | 2010-02-17 08:32:53 -0500 (Wed, 17 Feb 2010) | 1 line +Changed paths: + M /trunk/samtools/misc/samtools.pl + +Fixed filtering of SNPs near indels. Added min indel and SNP quality filter. +------------------------------------------------------------------------ +r535 | petulda | 2010-02-12 04:52:37 -0500 (Fri, 12 Feb 2010) | 1 line +Changed paths: + M /trunk/samtools/misc/sam2vcf.pl + +Print an error for pileups in simple format +------------------------------------------------------------------------ +r534 | lh3lh3 | 2010-02-11 14:01:41 -0500 (Thu, 11 Feb 2010) | 2 lines +Changed paths: + M /trunk/samtools/bam_plcmd.c + +added a hidden option in pileup to output the base position (for Erin) + +------------------------------------------------------------------------ +r533 | petulda | 2010-02-09 10:12:14 -0500 (Tue, 09 Feb 2010) | 1 line +Changed paths: + M /trunk/samtools/misc/sam2vcf.pl + +Added possibility to specify a custom column title for the data column +------------------------------------------------------------------------ +r532 | petulda | 2010-02-09 09:46:09 -0500 (Tue, 09 Feb 2010) | 1 line +Changed paths: + M /trunk/samtools/bam_plcmd.c + +Added the -d option to limit maximum depth for indels. +------------------------------------------------------------------------ +r531 | petulda | 2010-02-03 07:57:27 -0500 (Wed, 03 Feb 2010) | 1 line +Changed paths: + M /trunk/samtools/misc/sam2vcf.pl + +Added VCF header +------------------------------------------------------------------------ +r530 | lh3lh3 | 2010-02-01 09:13:19 -0500 (Mon, 01 Feb 2010) | 3 lines +Changed paths: + M /trunk/samtools/ChangeLog + M /trunk/samtools/bam_plcmd.c + M /trunk/samtools/bamtk.c + M /trunk/samtools/misc/samtools.pl + M /trunk/samtools/misc/wgsim.c + + * samtools-0.1.7-6 + * fixed a bug in faidx + ------------------------------------------------------------------------ r529 | jmarshall | 2010-01-11 18:51:49 -0500 (Mon, 11 Jan 2010) | 2 lines Changed paths: diff --git a/bam_aux.c b/bam_aux.c index 89e99f2..fbcd982 100644 --- a/bam_aux.c +++ b/bam_aux.c @@ -115,7 +115,7 @@ int bam_parse_region(bam_header_t *header, const char *str, int *ref_id, int *be *ref_id = kh_value(h, iter); if (i == k) { /* dump the whole sequence */ *begin = 0; *end = 1<<29; free(s); - return -1; + return 0; } for (p = s + i + 1; i != k; ++i) if (s[i] == '-') break; *begin = atoi(p); diff --git a/bam_index.c b/bam_index.c index 0d3e6fe..595ef04 100644 --- a/bam_index.c +++ b/bam_index.c @@ -475,6 +475,8 @@ int bam_index(int argc, char *argv[]) static inline int reg2bins(uint32_t beg, uint32_t end, uint16_t list[MAX_BIN]) { int i = 0, k; + if (beg >= end) return 0; + if (end >= 1u<<29) end = 1u<<29; --end; list[i++] = 0; for (k = 1 + (beg>>26); k <= 1 + (end>>26); ++k) list[i++] = k; @@ -496,7 +498,6 @@ struct __bam_iterf_t { int from_first; // read from the first record; no random access int tid, beg, end, n_off, i, finished; uint64_t curr_off; - const bam_index_t *idx; pair64_t *off; }; @@ -515,7 +516,7 @@ bam_iterf_t bam_iterf_query(const bam_index_t *idx, int tid, int beg, int end) if (end < beg) return 0; // initialize iter iter = calloc(1, sizeof(struct __bam_iterf_t)); - iter->idx = idx; iter->tid = tid, iter->beg = beg, iter->end = end; iter->i = -1; + iter->tid = tid, iter->beg = beg, iter->end = end; iter->i = -1; // bins = (uint16_t*)calloc(MAX_BIN, 2); n_bins = reg2bins(beg, end, bins); diff --git a/bam_pileup.c b/bam_pileup.c index 8a75829..d9a75da 100644 --- a/bam_pileup.c +++ b/bam_pileup.c @@ -232,6 +232,7 @@ const bam_pileup1_t *bam_plp_auto(bam_plp_t iter, int *_tid, int *_pos, int *_n_ if ((plp = bam_plp_next(iter, _tid, _pos, _n_plp)) != 0) return plp; else { *_n_plp = 0; + if (iter->is_eof) return 0; while (iter->func(iter->data, iter->b) >= 0) { if (bam_plp_push(iter, iter->b) < 0) { *_n_plp = -1; @@ -239,6 +240,8 @@ const bam_pileup1_t *bam_plp_auto(bam_plp_t iter, int *_tid, int *_pos, int *_n_ } if ((plp = bam_plp_next(iter, _tid, _pos, _n_plp)) != 0) return plp; } + bam_plp_push(iter, 0); + if ((plp = bam_plp_next(iter, _tid, _pos, _n_plp)) != 0) return plp; return 0; } } @@ -353,7 +356,7 @@ void bam_mplp_destroy(bam_mplp_t iter) { int i; for (i = 0; i < iter->n; ++i) bam_plp_destroy(iter->iter[i]); - free(iter->pos); free(iter->n_plp); free(iter->plp); + free(iter->iter); free(iter->pos); free(iter->n_plp); free(iter->plp); free(iter); } @@ -367,7 +370,7 @@ int bam_mplp_auto(bam_mplp_t iter, int *_tid, int *_pos, int *n_plp, const bam_p iter->plp[i] = bam_plp_auto(iter->iter[i], &tid, &pos, &iter->n_plp[i]); iter->pos[i] = (uint64_t)tid<<32 | pos; } - if (iter->pos[i] < new_min) new_min = iter->pos[i]; + if (iter->plp[i] && iter->pos[i] < new_min) new_min = iter->pos[i]; } iter->min = new_min; if (new_min == (uint64_t)-1) return 0; diff --git a/bam_plcmd.c b/bam_plcmd.c index 92023b8..d34cd26 100644 --- a/bam_plcmd.c +++ b/bam_plcmd.c @@ -163,18 +163,18 @@ static void pileup_seq(const bam_pileup1_t *p, int pos, int ref_len, const char if (p->is_head) printf("^%c", p->b->core.qual > 93? 126 : p->b->core.qual + 33); if (!p->is_del) { int j, rb, c = bam_nt16_rev_table[bam1_seqi(bam1_seq(p->b), p->qpos)]; - rb = (ref && (int)pos < ref_len)? ref[pos] : 'N'; + rb = (ref && pos < ref_len)? ref[pos] : 'N'; if (c == '=' || toupper(c) == toupper(rb)) c = bam1_strand(p->b)? ',' : '.'; else c = bam1_strand(p->b)? tolower(c) : toupper(c); putchar(c); if (p->indel > 0) { - putchar('+'); putw(p->indel, stdout); + printf("+%d", p->indel); for (j = 1; j <= p->indel; ++j) { c = bam_nt16_rev_table[bam1_seqi(bam1_seq(p->b), p->qpos + j)]; putchar(bam1_strand(p->b)? tolower(c) : toupper(c)); } } else if (p->indel < 0) { - putw(p->indel, stdout); + printf("%d", p->indel); for (j = 1; j <= -p->indel; ++j) { c = (ref && (int)pos+j < ref_len)? ref[pos+j] : 'N'; putchar(bam1_strand(p->b)? tolower(c) : toupper(c)); @@ -446,36 +446,81 @@ int bam_pileup(int argc, char *argv[]) * mpileup * ***********/ -static int mpileup(int n, char **fn) +typedef struct { + char *reg; + faidx_t *fai; +} mplp_conf_t; + +typedef struct { + bamFile fp; + bam_iterf_t iter; +} mplp_aux_t; + +static int mplp_func(void *data, bam1_t *b) { - bamFile *fp; - int i, tid, pos, *n_plp; + mplp_aux_t *ma = (mplp_aux_t*)data; + if (ma->iter) return bam_iterf_read(ma->fp, ma->iter, b); + return bam_read1(ma->fp, b); +} + +static int mpileup(mplp_conf_t *conf, int n, char **fn) +{ + mplp_aux_t **data; + int i, tid, pos, *n_plp, beg0 = 0, end0 = 1u<<29, ref_len, ref_tid; const bam_pileup1_t **plp; bam_mplp_t iter; bam_header_t *h = 0; - fp = calloc(n, sizeof(void*)); + char *ref; + // allocate + data = calloc(n, sizeof(void*)); plp = calloc(n, sizeof(void*)); n_plp = calloc(n, sizeof(int*)); + // read the header and initialize data for (i = 0; i < n; ++i) { bam_header_t *h_tmp; - fp[i] = bam_open(fn[i], "r"); - h_tmp = bam_header_read(fp[i]); + data[i] = calloc(1, sizeof(mplp_aux_t)); + data[i]->fp = bam_open(fn[i], "r"); + h_tmp = bam_header_read(data[i]->fp); + if (conf->reg) { + int beg, end; + bam_index_t *idx; + idx = bam_index_load(fn[i]); + if (idx == 0) { + fprintf(stderr, "[%s] fail to load index for %d-th input.\n", __func__, i+1); + exit(1); + } + if (bam_parse_region(h_tmp, conf->reg, &tid, &beg, &end) < 0) { + fprintf(stderr, "[%s] malformatted region or wrong seqname for %d-th input.\n", __func__, i+1); + exit(1); + } + if (i == 0) beg0 = beg, end0 = end; + data[i]->iter = bam_iterf_query(idx, tid, beg, end); + bam_index_destroy(idx); + } if (i == 0) h = h_tmp; else { // FIXME: to check consistency bam_header_destroy(h_tmp); } } - iter = bam_mplp_init(n, bam_read1, fp); + // mpileup + 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) { - printf("%s\t%d\tN", h->target_name[tid], pos + 1); + if (conf->reg && (pos < beg0 || pos >= end0)) continue; // out of the region requested + if (tid != ref_tid) { + free(ref); + 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, 0, 0); + 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; @@ -489,17 +534,35 @@ static int mpileup(int n, char **fn) } bam_mplp_destroy(iter); bam_header_destroy(h); - for (i = 0; i < n; ++i) bam_close(fp[i]); - free(fp); free(plp); + for (i = 0; i < n; ++i) { + bam_close(data[i]->fp); + if (data[i]->iter) bam_iterf_destroy(data[i]->iter); + free(data[i]); + } + free(data); free(plp); free(ref); free(n_plp); return 0; } 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:")) >= 0) { + switch (c) { + case 'f': + mplp.fai = fai_load(optarg); + if (mplp.fai == 0) return 1; + break; + case 'r': mplp.reg = strdup(optarg); + } + } if (argc == 1) { - fprintf(stderr, "Usage: samtools mpileup [ [...]]\n"); + fprintf(stderr, "Usage: samtools mpileup [-r reg] [-f in.fa] in1.bam [in2.bam [...]]\n"); return 1; } - mpileup(argc - 1, argv + 1); + mpileup(&mplp, argc - optind, argv + optind); + free(mplp.reg); + if (mplp.fai) fai_destroy(mplp.fai); return 0; } diff --git a/bamtk.c b/bamtk.c index 0594b11..ebd1b9e 100644 --- a/bamtk.c +++ b/bamtk.c @@ -9,7 +9,7 @@ #endif #ifndef PACKAGE_VERSION -#define PACKAGE_VERSION "0.1.7-14 (r591)" +#define PACKAGE_VERSION "0.1.7-15 (r592)" #endif int bam_taf2baf(int argc, char *argv[]); -- 2.39.2