examples/chk_indel.c

   1 /* To compile, copy this file to the samtools source code directory and compile with:
   2      gcc -g -O2 -Wall chk_indel_rg.c -o chk_indel_rg -Wall -I. -L. -lbam -lz
   3 */
   4
   5 #include <string.h>
   6 #include "bam.h"
   7
   8 typedef struct {
   9         long cnt[4]; // short:ins, short:del, long:ins, long:del
  10 } rgcnt_t;
  11
  12 #include "khash.h"
  13 KHASH_MAP_INIT_STR(rgcnt, rgcnt_t)
  14
  15 #define MAX_LEN 127
  16 #define Q_THRES 10
  17 #define L_THRES 6 // short: <=L_THRES; otherwise long
  18
  19 int main(int argc, char *argv[])
  20 {
  21         bamFile fp;
  22         bam1_t *b;
  23         int i, x;
  24         khash_t(rgcnt) *h;
  25         khint_t k;
  26
  27         if (argc == 1) {
  28                 fprintf(stderr, "Usage: chk_indel_rg <in.bam>\n\n");
  29                 fprintf(stderr, "Output: filename, RG, #ins-in-short-homopolymer, #del-in-short, #ins-in-long, #del-in-long\n");
  30                 return 1;
  31         }
  32
  33         h = kh_init(rgcnt);
  34         fp = bam_open(argv[1], "r");
  35         bam_header_destroy(bam_header_read(fp)); // we do not need the header
  36         b = bam_init1();
  37
  38         while (bam_read1(fp, b) >= 0) {
  39                 if (b->core.n_cigar >= 3 && b->core.qual >= Q_THRES) {
  40                         const uint8_t *seq;
  41                         const uint32_t *cigar = bam1_cigar(b);
  42                         char *rg;
  43                         for (i = 0; i < b->core.n_cigar; ++i) // check if there are 1bp indels
  44                                 if (bam_cigar_oplen(cigar[i]) == 1 && (bam_cigar_op(cigar[i]) == BAM_CDEL || bam_cigar_op(cigar[i]) == BAM_CINS))
  45                                         break;
  46                         if (i == b->core.n_cigar) continue; // no 1bp ins or del
  47                         if ((rg = (char*)bam_aux_get(b, "RG")) == 0) continue; // no RG tag
  48                         seq = bam1_seq(b);
  49                         for (i = x = 0; i < b->core.n_cigar; ++i) {
  50                                 int op = bam_cigar_op(cigar[i]);
  51                                 if (bam_cigar_oplen(cigar[i]) == 1 && (op == BAM_CDEL || op == BAM_CINS)) {
  52                                         int c, j, hrun, which;
  53                                         c = bam1_seqi(seq, x);
  54                                         for (j = x + 1, hrun = 0; j < b->core.l_qseq; ++j, ++hrun) // calculate the hompolymer run length
  55                                                 if (bam1_seqi(seq, j) != c) break;
  56                                         k = kh_get(rgcnt, h, rg + 1);
  57                                         if (k == kh_end(h)) { // absent
  58                                                 char *key = strdup(rg + 1);
  59                                                 k = kh_put(rgcnt, h, key, &c);
  60                                                 memset(&kh_val(h, k), 0, sizeof(rgcnt_t));
  61                                         }
  62                                         which = (hrun <= L_THRES? 0 : 1)<<1 | (op == BAM_CINS? 0 : 1);
  63                                         ++kh_val(h, k).cnt[which];
  64                                 }
  65                                 if (bam_cigar_type(op)&1) ++x;
  66                         }
  67                 }
  68         }
  69
  70         for (k = 0; k != kh_end(h); ++k) {
  71                 if (!kh_exist(h, k)) continue;
  72                 printf("%s\t%s", argv[1], kh_key(h, k));
  73                 for (i = 0; i < 4; ++i)
  74                         printf("\t%ld", kh_val(h, k).cnt[i]);
  75                 putchar('\n');
  76                 free((char*)kh_key(h, k));
  77         }
  78
  79         bam_destroy1(b);
  80         bam_close(fp);
  81         kh_destroy(rgcnt, h);
  82         return 0;
  83 }