samtools/sample.c.pysam.c

   1 #include "pysam.h"
   2
   3 #include <stdlib.h>
   4 #include <string.h>
   5 #include "sample.h"
   6 #include "khash.h"
   7 KHASH_MAP_INIT_STR(sm, int)
   8
   9 bam_sample_t *bam_smpl_init(void)
  10 {
  11         bam_sample_t *s;
  12         s = calloc(1, sizeof(bam_sample_t));
  13         s->rg2smid = kh_init(sm);
  14         s->sm2id = kh_init(sm);
  15         return s;
  16 }
  17
  18 void bam_smpl_destroy(bam_sample_t *sm)
  19 {
  20         int i;
  21         khint_t k;
  22         khash_t(sm) *rg2smid = (khash_t(sm)*)sm->rg2smid;
  23         if (sm == 0) return;
  24         for (i = 0; i < sm->n; ++i) free(sm->smpl[i]);
  25         free(sm->smpl);
  26         for (k = kh_begin(rg2smid); k != kh_end(rg2smid); ++k)
  27                 if (kh_exist(rg2smid, k)) free((char*)kh_key(rg2smid, k));
  28         kh_destroy(sm, sm->rg2smid);
  29         kh_destroy(sm, sm->sm2id);
  30         free(sm);
  31 }
  32
  33 static void add_pair(bam_sample_t *sm, khash_t(sm) *sm2id, const char *key, const char *val)
  34 {
  35         khint_t k_rg, k_sm;
  36         int ret;
  37         khash_t(sm) *rg2smid = (khash_t(sm)*)sm->rg2smid;
  38         k_rg = kh_get(sm, rg2smid, key);
  39         if (k_rg != kh_end(rg2smid)) return; // duplicated @RG-ID
  40         k_rg = kh_put(sm, rg2smid, strdup(key), &ret);
  41         k_sm = kh_get(sm, sm2id, val);
  42         if (k_sm == kh_end(sm2id)) { // absent
  43                 if (sm->n == sm->m) {
  44                         sm->m = sm->m? sm->m<<1 : 1;
  45                         sm->smpl = realloc(sm->smpl, sizeof(void*) * sm->m);
  46                 }
  47                 sm->smpl[sm->n] = strdup(val);
  48                 k_sm = kh_put(sm, sm2id, sm->smpl[sm->n], &ret);
  49                 kh_val(sm2id, k_sm) = sm->n++;
  50         }
  51         kh_val(rg2smid, k_rg) = kh_val(sm2id, k_sm);
  52 }
  53
  54 int bam_smpl_add(bam_sample_t *sm, const char *fn, const char *txt)
  55 {
  56         const char *p = txt, *q, *r;
  57         kstring_t buf, first_sm;
  58         int n = 0;
  59         khash_t(sm) *sm2id = (khash_t(sm)*)sm->sm2id;
  60         if (txt == 0) {
  61                 add_pair(sm, sm2id, fn, fn);
  62                 return 0;
  63         }
  64         memset(&buf, 0, sizeof(kstring_t));
  65     memset(&first_sm, 0, sizeof(kstring_t));
  66         while ((q = strstr(p, "@RG")) != 0) {
  67                 p = q + 3;
  68                 r = q = 0;
  69                 if ((q = strstr(p, "\tID:")) != 0) q += 4;
  70                 if ((r = strstr(p, "\tSM:")) != 0) r += 4;
  71                 if (r && q) {
  72                         char *u, *v;
  73                         int oq, or;
  74                         for (u = (char*)q; *u && *u != '\t' && *u != '\n'; ++u);
  75                         for (v = (char*)r; *v && *v != '\t' && *v != '\n'; ++v);
  76                         oq = *u; or = *v; *u = *v = '\0';
  77                         buf.l = 0; kputs(fn, &buf); kputc('/', &buf); kputs(q, &buf);
  78                         add_pair(sm, sm2id, buf.s, r);
  79             if ( !first_sm.s )
  80                 kputs(r,&first_sm);
  81                         *u = oq; *v = or;
  82                 } else break;
  83                 p = q > r? q : r;
  84                 ++n;
  85         }
  86         if (n == 0) add_pair(sm, sm2id, fn, fn);
  87     // If there is only one RG tag present in the header and reads are not annotated, don't refuse to work but
  88     //  use the tag instead.
  89     else if ( n==1 && first_sm.s )
  90         add_pair(sm,sm2id,fn,first_sm.s);
  91     if ( first_sm.s )
  92         free(first_sm.s);
  93
  94 //      add_pair(sm, sm2id, fn, fn);
  95         free(buf.s);
  96         return 0;
  97 }
  98
  99 int bam_smpl_rg2smid(const bam_sample_t *sm, const char *fn, const char *rg, kstring_t *str)
 100 {
 101         khint_t k;
 102         khash_t(sm) *rg2smid = (khash_t(sm)*)sm->rg2smid;
 103         if (rg) {
 104                 str->l = 0;
 105                 kputs(fn, str); kputc('/', str); kputs(rg, str);
 106                 k = kh_get(sm, rg2smid, str->s);
 107         } else k = kh_get(sm, rg2smid, fn);
 108         return k == kh_end(rg2smid)? -1 : kh_val(rg2smid, k);
 109 }