* samtools-0.1.4-8 (r338)

 * parse the @RG header lines and allow to choose library at the "samtools view"
   command line

diff --git a/bam.c b/bam.c
index 22e09cf15c97a9035557926a1eda4c8e07792586..e3cb6bd0526ff7e97403cf0b38882bda51075a7a 100644 (file)
--- a/bam.c
+++ b/bam.c
@@ -79,6 +79,7 @@ void bam_header_destroy(bam_header_t *header)
        }
        free(header->text);
 #ifndef BAM_NO_HASH
+       if (header->rg2lib) bam_strmap_destroy(header->rg2lib);
        bam_destroy_header_hash(header);
 #endif
        free(header);

diff --git a/bam.h b/bam.h
index 9330075ed10bcbcd67e3fba245f87909b56b42ef..c635e72767e31852ddf8b2759a3d9493d13e58bc 100644 (file)
--- a/bam.h
+++ b/bam.h
@@ -100,7 +100,7 @@ typedef struct {
        int32_t n_targets;
        char **target_name;
        uint32_t *target_len;
-       void *hash;
+       void *hash, *rg2lib;
        int l_text;
        char *text;
 } bam_header_t;
@@ -317,9 +317,16 @@ extern "C" {
 
        bam_header_t *sam_header_read(tamFile fp);
        int sam_header_parse(bam_header_t *h);
+       int sam_header_parse_rg(bam_header_t *h);
 
 #define sam_write1(header, b) bam_view1(header, b)
 
+       int bam_strmap_put(void *strmap, const char *rg, const char *lib);
+       const char *bam_strmap_get(const void *strmap, const char *rg);
+       void *bam_strmap_dup(const void*);
+       void *bam_strmap_init();
+       void bam_strmap_destroy(void *strmap);
+
        /*!
          @abstract Initialize a header structure.
          @return   the pointer to the header structure
@@ -593,7 +600,7 @@ extern "C" {
 
        void bam_aux_append(bam1_t *b, const char tag[2], char type, int len, uint8_t *data);
        // uint8_t *bam_aux_get_core(bam1_t *b, const char tag[2]); // an alias of bam_aux_get()
-       uint8_t *bam_aux_get(bam1_t *b, const char tag[2]);
+       uint8_t *bam_aux_get(const bam1_t *b, const char tag[2]);
        int32_t bam_aux2i(const uint8_t *s);
        float bam_aux2f(const uint8_t *s);
        double bam_aux2d(const uint8_t *s);

diff --git a/bam_aux.c b/bam_aux.c
index 41be81910d89478dc160fada8fe2b4c38f8052c5..137f9afffdd016347a48cd85df7a6c0cfd2cebf0 100644 (file)
--- a/bam_aux.c
+++ b/bam_aux.c
@@ -1,7 +1,9 @@
 #include <ctype.h>
 #include "bam.h"
 #include "khash.h"
+typedef char *str_p;
 KHASH_MAP_INIT_STR(s, int)
+KHASH_MAP_INIT_STR(r2l, str_p)
 
 void bam_aux_append(bam1_t *b, const char tag[2], char type, int len, uint8_t *data)
 {
@@ -23,7 +25,7 @@ uint8_t *bam_aux_get_core(bam1_t *b, const char tag[2])
        return bam_aux_get(b, tag);
 }
 */
-uint8_t *bam_aux_get(bam1_t *b, const char tag[2])
+uint8_t *bam_aux_get(const bam1_t *b, const char tag[2])
 {
        uint8_t *s;
        int y = tag[0]<<8 | tag[1];
@@ -158,6 +160,70 @@ char *bam_aux2Z(const uint8_t *s)
        else return 0;
 }
 
+/******************
+ * rg2lib related *
+ ******************/
+
+int bam_strmap_put(void *rg2lib, const char *rg, const char *lib)
+{
+       int ret;
+       khint_t k;
+       khash_t(r2l) *h = (khash_t(r2l)*)rg2lib;
+       char *key = strdup(rg);
+       k = kh_put(r2l, h, key, &ret);
+       if (ret) kh_val(h, k) = strdup(lib);
+       else {
+               fprintf(stderr, "[bam_rg2lib_put] duplicated @RG ID: %s\n", rg);
+               free(key);
+       }
+       return 0;
+}
+
+const char *bam_strmap_get(const void *rg2lib, const char *rg)
+{
+       const khash_t(r2l) *h = (const khash_t(r2l)*)rg2lib;
+       khint_t k;
+       k = kh_get(r2l, h, rg);
+       if (k != kh_end(h)) return (const char*)kh_val(h, k);
+       else return 0;
+}
+
+void *bam_strmap_dup(const void *rg2lib)
+{
+       const khash_t(r2l) *h = (const khash_t(r2l)*)rg2lib;
+       khash_t(r2l) *g = kh_init(r2l);
+       khint_t k, l;
+       int ret;
+       for (k = kh_begin(h); k < kh_end(h); ++k) {
+               if (kh_exist(h, k)) {
+                       char *key = strdup(kh_key(h, k));
+                       l = kh_put(r2l, g, key, &ret);
+                       kh_val(g, l) = strdup(kh_val(h, k));
+               }
+       }
+       return g;
+}
+
+void *bam_strmap_init()
+{
+       return (void*)kh_init(r2l);
+}
+
+void bam_strmap_destroy(void *rg2lib)
+{
+       khash_t(r2l) *h = (khash_t(r2l)*)rg2lib;
+       khint_t k;
+       if (h == 0) return;
+       for (k = kh_begin(h); k < kh_end(h); ++k) {
+               if (kh_exist(h, k)) {
+                       free((char*)kh_key(h, k)); free(kh_val(h, k));
+               }
+       }
+       kh_destroy(r2l, h);
+}
+
+/*** The following routines were implemented by Nils Homer for color-space support in tview ***/
+
 char bam_aux_getCSi(bam1_t *b, int i)
 {
        uint8_t *c = bam_aux_get(b, "CS");

diff --git a/bam_import.c b/bam_import.c
index 77e324319525b529ac456aad2e9d9f84769e7b2d..c08547345056a356b6b07b50e931c9baff05f48c 100644 (file)
--- a/bam_import.c
+++ b/bam_import.c
@@ -5,6 +5,7 @@
 #include <stdlib.h>
 #include <unistd.h>
 #include <assert.h>
+#include "kstring.h"
 #include "bam.h"
 #include "kseq.h"
 #include "khash.h"
@@ -146,6 +147,55 @@ static inline void append_text(bam_header_t *header, kstring_t *str)
        header->text[header->l_text] = 0;
 }
 
+int sam_header_parse_rg(bam_header_t *h)
+{
+       kstring_t *rgid, *rglib;
+       char *p, *q, *s, *r;
+       int n = 0;
+
+       bam_strmap_destroy(h->rg2lib);
+       // parse @RG lines
+       h->rg2lib = bam_strmap_init();
+       rgid = calloc(1, sizeof(kstring_t));
+       rglib = calloc(1, sizeof(kstring_t));
+       s = h->text;
+       while ((s = strstr(s, "@RG")) != 0) {
+               if (rgid->l && rglib->l) {
+                       bam_strmap_put(h->rg2lib, rgid->s, rglib->s);
+                       ++n;
+               }
+               rgid->l = rglib->l = 0;
+               s += 3;
+               r = s;
+               if ((p = strstr(s, "ID:")) != 0) {
+                       q = p + 3;
+                       for (p = q; *p && *p != '\t' && *p != '\r' && *p != '\n'; ++p);
+                       kputsn(q, p - q, rgid);
+               } else {
+                       fprintf(stderr, "[bam_header_parse] missing ID tag in @RG lines.\n");
+                       break;
+               }
+               if (r < p) r = p;
+               if ((p = strstr(s, "LB:")) != 0) {
+                       q = p + 3;
+                       for (p = q; *p && *p != '\t' && *p != '\r' && *p != '\n'; ++p);
+                       kputsn(q, p - q, rglib);
+               } else {
+                       fprintf(stderr, "[bam_header_parse] missing LB tag in @RG lines.\n");
+                       break;
+               }
+               if (r < p) r = p;
+               s = r + 3;
+       }
+       if (rgid->l && rglib->l) {
+               bam_strmap_put(h->rg2lib, rgid->s, rglib->s);
+               ++n;
+       }
+       free(rgid->s); free(rgid);
+       free(rglib->s); free(rglib);
+       return n;
+}
+
 int sam_header_parse(bam_header_t *h)
 {
        int i;
@@ -183,6 +233,7 @@ int sam_header_parse(bam_header_t *h)
                s = r + 3;
                ++i;
        }
+       sam_header_parse_rg(h);
        return h->n_targets;
 
 header_err_ret:

diff --git a/bamtk.c b/bamtk.c
index 878562cca946bb417779456d329230eb2c1d9da4..852213e3580c90de58052537ed19f003a0340de8 100644 (file)
--- a/bamtk.c
+++ b/bamtk.c
@@ -3,7 +3,7 @@
 #include "bam.h"
 
 #ifndef PACKAGE_VERSION
-#define PACKAGE_VERSION "0.1.4-7 (r337)"
+#define PACKAGE_VERSION "0.1.4-8 (r338)"
 #endif
 
 int bam_taf2baf(int argc, char *argv[]);

diff --git a/kstring.h b/kstring.h
index 47d63f7690d7e146502856050b38a3699b810ba0..221ade2472655ba842477b1249feb7ca7209631c 100644 (file)
--- a/kstring.h
+++ b/kstring.h
@@ -19,19 +19,24 @@ typedef struct __kstring_t {
 int ksprintf(kstring_t *s, const char *fmt, ...);
 int ksplit_core(char *s, int delimiter, int *_max, int **_offsets);
 
-static inline int kputs(const char *p, kstring_t *s)
+static inline int kputsn(const char *p, int l, kstring_t *s)
 {
-       int l = strlen(p);
        if (s->l + l + 1 >= s->m) {
                s->m = s->l + l + 2;
                kroundup32(s->m);
                s->s = (char*)realloc(s->s, s->m);
        }
-       strcpy(s->s + s->l, p);
+       strncpy(s->s + s->l, p, l);
        s->l += l;
+       s->s[s->l] = 0;
        return l;
 }
 
+static inline int kputs(const char *p, kstring_t *s)
+{
+       return kputsn(p, strlen(p), s);
+}
+
 static inline int kputc(int c, kstring_t *s)
 {
        if (s->l + 1 >= s->m) {

diff --git a/sam.c b/sam.c
index 04be5c57745a8036dce8b2e3ede563f38aa8cdbb..4c16b02fb87de19426f50dc3e1573b3f58801a26 100644 (file)
--- a/sam.c
+++ b/sam.c
@@ -19,6 +19,7 @@ bam_header_t *bam_header_dup(const bam_header_t *h0)
                h->target_len[i] = h0->target_len[i];
                h->target_name[i] = strdup(h0->target_name[i]);
        }
+       if (h0->rg2lib) h->rg2lib = bam_strmap_dup(h0->rg2lib);
        return h;
 }
 
@@ -46,6 +47,7 @@ samfile_t *samopen(const char *fn, const char *mode, const void *aux)
                                        fprintf(stderr, "[samopen] no @SQ lines in the header.\n");
                        } else fprintf(stderr, "[samopen] SAM header is present: %d sequences.\n", fp->header->n_targets);
                }
+               sam_header_parse_rg(fp->header);
        } else if (mode[0] == 'w') { // write
                fp->header = bam_header_dup((const bam_header_t*)aux);
                if (mode[1] == 'b') { // binary

diff --git a/sam_view.c b/sam_view.c
index 48586106aabe3fb9bdad2f7c48d7acba8e3f36ea..50192f17f647bf64dc1b060343095d37ce120218 100644 (file)
--- a/sam_view.c
+++ b/sam_view.c
@@ -5,14 +5,26 @@
 #include "sam.h"
 
 static int g_min_mapQ = 0, g_flag_on = 0, g_flag_off = 0;
+static char *g_library;
 
-#define __g_skip_aln(b) (((b)->core.qual < g_min_mapQ) || ((b->core.flag & g_flag_on) != g_flag_on) \
-                                                || (b->core.flag & g_flag_off))
+static inline int __g_skip_aln(const bam_header_t *h, const bam1_t *b)
+{
+       if (b->core.qual < g_min_mapQ || ((b->core.flag & g_flag_on) != g_flag_on) || (b->core.flag & g_flag_off))
+               return 1;
+       if (g_library) {
+               uint8_t *s = bam_aux_get(b, "RG");
+               if (s) {
+                       const char *p = bam_strmap_get(h->rg2lib, s + 1);
+                       return (p && strcmp(p, g_library) == 0)? 0 : 1;
+               } else return 1;
+       } else return 0;
+}
 
 // callback function for bam_fetch()
 static int view_func(const bam1_t *b, void *data)
 {
-       if (!__g_skip_aln(b)) samwrite((samfile_t*)data, b);
+       if (!__g_skip_aln(((samfile_t*)data)->header, b))
+               samwrite((samfile_t*)data, b);
        return 0;
 }
 
@@ -26,7 +38,7 @@ int main_samview(int argc, char *argv[])
 
        /* parse command-line options */
        strcpy(in_mode, "r"); strcpy(out_mode, "w");
-       while ((c = getopt(argc, argv, "Sbt:hHo:q:f:F:u")) >= 0) {
+       while ((c = getopt(argc, argv, "Sbt:hHo:q:f:F:ul:")) >= 0) {
                switch (c) {
                case 'S': is_bamin = 0; break;
                case 'b': is_bamout = 1; break;
@@ -38,6 +50,7 @@ int main_samview(int argc, char *argv[])
                case 'F': g_flag_off = strtol(optarg, 0, 0); break;
                case 'q': g_min_mapQ = atoi(optarg); break;
                case 'u': is_uncompressed = 1; break;
+               case 'l': g_library = strdup(optarg); break;
                default: return usage();
                }
        }
@@ -64,7 +77,7 @@ int main_samview(int argc, char *argv[])
                bam1_t *b = bam_init1();
                int r;
                while ((r = samread(in, b)) >= 0) // read one alignment from `in'
-                       if (!__g_skip_aln(b))
+                       if (!__g_skip_aln(in->header, b))
                                samwrite(out, b); // write the alignment to `out'
                if (r < -1) fprintf(stderr, "[main_samview] truncated file.\n");
                bam_destroy1(b);
@@ -91,7 +104,7 @@ int main_samview(int argc, char *argv[])
 
 view_end:
        // close files, free and return
-       free(fn_list); free(fn_out);
+       free(fn_list); free(fn_out); free(g_library);
        samclose(in);
        samclose(out);
        return ret;