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