bam_calnh.c

   1 #include <stdio.h>
   2 #include <stdlib.h>
   3
   4 #include "sam.h"
   5 #include "khash.h"
   6
   7 KHASH_MAP_INIT_STR(QNH, int)
   8
   9 int
  10 bam_calnh_core(samfile_t *in, samfile_t *out)
  11 {
  12   bam1_t *b;
  13   khint_t k;
  14   uint8_t *oldnh;
  15   char *qn = NULL;
  16   khash_t(QNH) *qnh;
  17   int ret;
  18   int64_t savepos;
  19   uint64_t val;
  20
  21   savepos = bam_tell(in->x.bam);
  22
  23   if (savepos < 0)
  24   {
  25     fprintf(stderr, "[bam_calnh_core] bam_tell failure\n");
  26     return 1;
  27   }
  28
  29   qnh = kh_init(QNH);
  30   b = bam_init1();
  31
  32   while (samread(in, b) >= 0)
  33   {
  34     qn = strdup(bam1_qname(b));
  35     k = kh_put(QNH, qnh, qn, &ret);
  36
  37     if (ret == 0)
  38     {
  39       kh_val(qnh, k) += 1;
  40       free(qn);
  41     }
  42
  43     else
  44       kh_val(qnh, k) = 1;
  45   }
  46
  47   if ((bam_seek(in->x.bam, savepos, SEEK_SET)) < 0)
  48   {
  49     fprintf(stderr, "[bam_calnh_core] bam_seek failure\n");
  50     return 1;
  51   }
  52
  53
  54   while (samread(in, b) >= 0)
  55   {
  56     k = kh_get(QNH, qnh, bam1_qname(b));
  57     oldnh = bam_aux_get(b, "NH");
  58
  59     if (!(oldnh))
  60     {
  61       val = kh_val(qnh, k);
  62       bam_aux_append(b, "NH", 'i', 4, (uint8_t *)&val);
  63     }
  64
  65     samwrite(out, b);
  66   }
  67
  68   for (k = kh_begin(qnh); k != kh_end(qnh); ++k)
  69     if (kh_exist(qnh, k))
  70       free((char *)kh_key(qnh, k));
  71
  72   kh_destroy(QNH, qnh);
  73   bam_destroy1(b);
  74   return 0;
  75 }
  76
  77 void
  78 usage(char **argv)
  79 {
  80 #ifdef _STANDALONE
  81     fprintf(stderr, "\nUsage:  %s <input.bam> <output.bam>\n", argv[0]);
  82 #else
  83     fprintf(stderr, "\nUsage:  samtools calnh <input.bam> <output.bam>\n");
  84 #endif /* _STANDALONE */
  85     fprintf(stderr, "    specify input.bam or output.bam as '-' to use pipes\n\n");
  86 }
  87
  88 int
  89 bam_calnh(int argc, char **argv)
  90 {
  91   samfile_t *in;
  92   samfile_t *out;
  93   int ret;
  94
  95   if (argc < 3)
  96   {
  97     usage(argv);
  98     return 1;
  99   }
 100
 101   in = samopen(argv[1], "rb", 0);
 102   out = samopen(argv[2], "wb", in->header);
 103
 104   if ((in == 0) || (out == 0))
 105   {
 106     fprintf(stderr, "[bam_calnh] fail to read/write input files\n");
 107     return 1;
 108   }
 109
 110   ret = bam_calnh_core(in, out);
 111   samclose(in);
 112   samclose(out);
 113   return ret;
 114 }
 115
 116 #ifdef _STANDALONE
 117 int
 118 main(int argc, char **argv)
 119 {
 120   return bam_calnh(argc, argv);
 121 }
 122 #endif /* _STANDALONE */