bam_stat.c

   1 #include <unistd.h>
   2 #include "bam.h"
   3
   4 typedef struct {
   5         long long n_reads, n_mapped, n_pair_all, n_pair_map, n_pair_good;
   6         long long n_sgltn, n_read1, n_read2;
   7         long long n_qcfail, n_dup;
   8         long long n_diffchr, n_diffhigh;
   9 } bam_flagstat_t;
  10
  11 #define flagstat_loop(s, c) do {                                                                                \
  12                 ++(s)->n_reads;                                                                                                 \
  13                 if ((c)->flag & BAM_FPAIRED) {                                                                  \
  14                         ++(s)->n_pair_all;                                                                                      \
  15                         if ((c)->flag & BAM_FPROPER_PAIR) ++(s)->n_pair_good;           \
  16                         if ((c)->flag & BAM_FREAD1) ++(s)->n_read1;                                     \
  17                         if ((c)->flag & BAM_FREAD2) ++(s)->n_read2;                                     \
  18                         if ((c)->flag & BAM_FMUNMAP) ++(s)->n_sgltn;                            \
  19                         if (!((c)->flag & BAM_FUNMAP) && !((c)->flag & BAM_FMUNMAP)) { \
  20                                 ++(s)->n_pair_map;                                                                              \
  21                                 if ((c)->mtid != (c)->tid) {                                                    \
  22                                         ++(s)->n_diffchr;                                                                       \
  23                                         if ((c)->qual >= 5) ++(s)->n_diffhigh;                          \
  24                                 }                                                                                                               \
  25                         }                                                                                                                       \
  26                 }                                                                                                                               \
  27                 if (!((c)->flag & BAM_FUNMAP)) ++(s)->n_mapped;                                 \
  28                 if ((c)->flag & BAM_FQCFAIL) ++(s)->n_qcfail;                                   \
  29                 if ((c)->flag & BAM_FDUP) ++(s)->n_dup;                                                 \
  30         } while (0)
  31
  32 bam_flagstat_t *bam_flagstat_core(bamFile fp)
  33 {
  34         bam_flagstat_t *s;
  35         bam1_t *b;
  36         bam1_core_t *c;
  37         int ret;
  38         s = (bam_flagstat_t*)calloc(1, sizeof(bam_flagstat_t));
  39         b = bam_init1();
  40         c = &b->core;
  41         while ((ret = bam_read1(fp, b)) >= 0)
  42                 flagstat_loop(s, c);
  43         bam_destroy1(b);
  44         if (ret != -1)
  45                 fprintf(stderr, "[bam_flagstat_core] Truncated file? Continue anyway.\n");
  46         return s;
  47 }
  48 int bam_flagstat(int argc, char *argv[])
  49 {
  50         bamFile fp;
  51         bam_header_t *header;
  52         bam_flagstat_t *s;
  53         if (argc == optind) {
  54                 fprintf(stderr, "Usage: samtools flagstat <in.bam>\n");
  55                 return 1;
  56         }
  57         fp = strcmp(argv[optind], "-")? bam_open(argv[optind], "r") : bam_dopen(fileno(stdin), "r");
  58         assert(fp);
  59         header = bam_header_read(fp);
  60         s = bam_flagstat_core(fp);
  61         printf("%lld in total\n", s->n_reads);
  62         printf("%lld QC failure\n", s->n_qcfail);
  63         printf("%lld duplicates\n", s->n_dup);
  64         printf("%lld mapped (%.2f%%)\n", s->n_mapped, (float)s->n_mapped / s->n_reads * 100.0);
  65         printf("%lld paired in sequencing\n", s->n_pair_all);
  66         printf("%lld read1\n", s->n_read1);
  67         printf("%lld read2\n", s->n_read2);
  68         printf("%lld properly paired (%.2f%%)\n", s->n_pair_good, (float)s->n_pair_good / s->n_pair_all * 100.0);
  69         printf("%lld with itself and mate mapped\n", s->n_pair_map);
  70         printf("%lld singletons (%.2f%%)\n", s->n_sgltn, (float)s->n_sgltn / s->n_pair_all * 100.0);
  71         printf("%lld with mate mapped to a different chr\n", s->n_diffchr);
  72         printf("%lld with mate mapped to a different chr (mapQ>=5)\n", s->n_diffhigh);
  73         free(s);
  74         bam_header_destroy(header);
  75         bam_close(fp);
  76         return 0;
  77 }