Mercurial > repos > lsong10 > psiclass

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

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author lsong10
date Fri, 26 Mar 2021 16:52:45 +0000
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:903fc43d6227
1 #include <zlib.h>
2 #include <stdio.h>
3 #include <ctype.h>
4 #include <stdlib.h>
5 #include <string.h>
6 #include <unistd.h>
7 #include "kstring.h"
8 #include "bgzf.h"
9 #include "bam.h"
10
11 #include "kseq.h"
12 KSTREAM_INIT(gzFile, gzread, 16384)
13
14 typedef struct {
15 bamFile fp;
16 bam_iter_t iter;
17 int min_mapQ;
18 } aux_t;
19
20 static int read_bam(void *data, bam1_t *b)
21 {
22 aux_t *aux = (aux_t*)data;
23 int ret = bam_iter_read(aux->fp, aux->iter, b);
24 if ((int)b->core.qual < aux->min_mapQ) b->core.flag |= BAM_FUNMAP;
25 return ret;
26 }
27
28 int main_bedcov(int argc, char *argv[])
29 {
30 extern void bam_init_header_hash(bam_header_t*);
31 gzFile fp;
32 kstring_t str;
33 kstream_t *ks;
34 bam_index_t **idx;
35 bam_header_t *h = 0;
36 aux_t **aux;
37 int *n_plp, dret, i, n, c, min_mapQ = 0;
38 int64_t *cnt;
39 const bam_pileup1_t **plp;
40
41 while ((c = getopt(argc, argv, "Q:")) >= 0) {
42 switch (c) {
43 case 'Q': min_mapQ = atoi(optarg); break;
44 }
45 }
46 if (optind + 2 > argc) {
47 fprintf(stderr, "Usage: samtools bedcov <in.bed> <in1.bam> [...]\n");
48 return 1;
49 }
50 memset(&str, 0, sizeof(kstring_t));
51 n = argc - optind - 1;
52 aux = calloc(n, sizeof(void*));
53 idx = calloc(n, sizeof(void*));
54 for (i = 0; i < n; ++i) {
55 aux[i] = calloc(1, sizeof(aux_t));
56 aux[i]->min_mapQ = min_mapQ;
57 aux[i]->fp = bam_open(argv[i+optind+1], "r");
58 idx[i] = bam_index_load(argv[i+optind+1]);
59 if (aux[i]->fp == 0 || idx[i] == 0) {
60 fprintf(stderr, "ERROR: fail to open index BAM file '%s'\n", argv[i+optind+1]);
61 return 2;
62 }
63 bgzf_set_cache_size(aux[i]->fp, 20);
64 if (i == 0) h = bam_header_read(aux[0]->fp);
65 }
66 bam_init_header_hash(h);
67 cnt = calloc(n, 8);
68
69 fp = gzopen(argv[optind], "rb");
70 ks = ks_init(fp);
71 n_plp = calloc(n, sizeof(int));
72 plp = calloc(n, sizeof(void*));
73 while (ks_getuntil(ks, KS_SEP_LINE, &str, &dret) >= 0) {
74 char *p, *q;
75 int tid, beg, end, pos;
76 bam_mplp_t mplp;
77
78 for (p = q = str.s; *p && *p != '\t'; ++p);
79 if (*p != '\t') goto bed_error;
80 *p = 0; tid = bam_get_tid(h, q); *p = '\t';
81 if (tid < 0) goto bed_error;
82 for (q = p = p + 1; isdigit(*p); ++p);
83 if (*p != '\t') goto bed_error;
84 *p = 0; beg = atoi(q); *p = '\t';
85 for (q = p = p + 1; isdigit(*p); ++p);
86 if (*p == '\t' || *p == 0) {
87 int c = *p;
88 *p = 0; end = atoi(q); *p = c;
89 } else goto bed_error;
90
91 for (i = 0; i < n; ++i) {
92 if (aux[i]->iter) bam_iter_destroy(aux[i]->iter);
93 aux[i]->iter = bam_iter_query(idx[i], tid, beg, end);
94 }
95 mplp = bam_mplp_init(n, read_bam, (void**)aux);
96 bam_mplp_set_maxcnt(mplp, 64000);
97 memset(cnt, 0, 8 * n);
98 while (bam_mplp_auto(mplp, &tid, &pos, n_plp, plp) > 0)
99 if (pos >= beg && pos < end)
100 for (i = 0; i < n; ++i) cnt[i] += n_plp[i];
101 for (i = 0; i < n; ++i) {
102 kputc('\t', &str);
103 kputl(cnt[i], &str);
104 }
105 puts(str.s);
106 bam_mplp_destroy(mplp);
107 continue;
108
109 bed_error:
110 fprintf(stderr, "Errors in BED line '%s'\n", str.s);
111 }
112 free(n_plp); free(plp);
113 ks_destroy(ks);
114 gzclose(fp);
115
116 free(cnt);
117 for (i = 0; i < n; ++i) {
118 if (aux[i]->iter) bam_iter_destroy(aux[i]->iter);
119 bam_index_destroy(idx[i]);
120 bam_close(aux[i]->fp);
121 free(aux[i]);
122 }
123 bam_header_destroy(h);
124 free(aux); free(idx);
125 free(str.s);
126 return 0;
127 }