]> git.donarmstrong.com Git - samtools.git/commitdiff
* samtools-0.1.7-15 (r592)
authorHeng Li <lh3@live.co.uk>
Sat, 12 Jun 2010 22:02:38 +0000 (22:02 +0000)
committerHeng Li <lh3@live.co.uk>
Sat, 12 Jun 2010 22:02:38 +0000 (22:02 +0000)
 * fixed a few minor memory leaks in the new pileup code
 * improved the functionality of mpileup

ChangeLog
bam_aux.c
bam_index.c
bam_pileup.c
bam_plcmd.c
bamtk.c

index 64df707683cb8c2ee3fb4a6cf7f0612482a37413..2e6f5354e1856a47094ab14181e6ae103291eb0b 100644 (file)
--- 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:
index 89e99f281adb652e144eea0a77074a1c6cec023a..fbcd9822b233dab64f6c861fb332846acb951ce7 100644 (file)
--- 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);
index 0d3e6fe8aca611b73b83bcb3fb7bdafd087f6753..595ef047733dc62282732c0dd6ece2d29537887f 100644 (file)
@@ -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);
index 8a75829cf6e92889ed72e2b23f08056abedb516d..d9a75da013eecf772336bf5a03ef9baa4851b41d 100644 (file)
@@ -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;
index 92023b87f92e4c293c7bb4409cac63373668c017..d34cd267a9edd363a32605c083fefe5ba512e8bb 100644 (file)
@@ -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 <in1.bam> [<in2.bam> [...]]\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 0594b119ba13ea8960fec6601c88ebf7fcf5128e..ebd1b9e2dc9aa55b8b3d45c9507861912d321274 100644 (file)
--- 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[]);