annotate PsiCLASS-1.0.2/samtools-0.1.19/bam_stat.c @ 0:903fc43d6227 draft default tip

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