Mercurial > repos > youngkim > ezbamqc

comparison ezBAMQC/src/htslib/test/test_view.c @ 0:dfa3745e5fd8

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author youngkim
date Thu, 24 Mar 2016 17:12:52 -0400
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:dfa3745e5fd8
1 /* test/test_view.c -- simple view tool, purely for use in a test harness.
2
3 Copyright (C) 2012 Broad Institute.
4 Copyright (C) 2013-2014 Genome Research Ltd.
5
6 Author: Heng Li <lh3@sanger.ac.uk>
7
8 Permission is hereby granted, free of charge, to any person obtaining a copy
9 of this software and associated documentation files (the "Software"), to deal
10 in the Software without restriction, including without limitation the rights
11 to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
12 copies of the Software, and to permit persons to whom the Software is
13 furnished to do so, subject to the following conditions:
14
15 The above copyright notice and this permission notice shall be included in
16 all copies or substantial portions of the Software.
17
18 THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
19 IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
20 FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
21 THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
22 LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
23 FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
24 DEALINGS IN THE SOFTWARE. */
25
26 #include <stdio.h>
27 #include <unistd.h>
28 #include <stdlib.h>
29 #include <string.h>
30
31 #include "cram/cram.h"
32
33 #include "htslib/sam.h"
34
35 typedef struct hts_opt {
36 enum cram_option opt;
37 union {
38 int i;
39 char *s;
40 } val;
41 struct hts_opt *next;
42 } hts_opt;
43
44 /*
45 * Parses arg and appends it to the option list.
46 * Returns 0 on success;
47 * -1 on failure.
48 */
49 int add_option(hts_opt **opts, char *arg) {
50 hts_opt *o, *t;
51 char *cp;
52
53 if (!(cp = strchr(arg, '=')))
54 cp = "1"; // assume boolean
55 else
56 *cp++ = 0;
57
58 if (!(o = malloc(sizeof(*o))))
59 return -1;
60
61 if (strcmp(arg, "DECODE_MD") == 0)
62 o->opt = CRAM_OPT_DECODE_MD, o->val.i = atoi(cp);
63 else if (strcmp(arg, "VERBOSITY") == 0)
64 o->opt = CRAM_OPT_VERBOSITY, o->val.i = atoi(cp);
65 else if (strcmp(arg, "SEQS_PER_SLICE") == 0)
66 o->opt = CRAM_OPT_SEQS_PER_SLICE, o->val.i = atoi(cp);
67 else if (strcmp(arg, "SLICES_PER_CONTAINER") == 0)
68 o->opt = CRAM_OPT_SLICES_PER_CONTAINER, o->val.i = atoi(cp);
69 else if (strcmp(arg, "EMBED_REF") == 0)
70 o->opt = CRAM_OPT_EMBED_REF, o->val.i = atoi(cp);
71 else if (strcmp(arg, "NO_REF") == 0)
72 o->opt = CRAM_OPT_NO_REF, o->val.i = atoi(cp);
73 else if (strcmp(arg, "IGNORE_MD5") == 0)
74 o->opt = CRAM_OPT_IGNORE_MD5, o->val.i = atoi(cp);
75 else if (strcmp(arg, "USE_BZIP2") == 0)
76 o->opt = CRAM_OPT_USE_BZIP2, o->val.i = atoi(cp);
77 else if (strcmp(arg, "USE_RANS") == 0)
78 o->opt = CRAM_OPT_USE_RANS, o->val.i = atoi(cp);
79 else if (strcmp(arg, "USE_LZMA") == 0)
80 o->opt = CRAM_OPT_USE_LZMA, o->val.i = atoi(cp);
81 else if (strcmp(arg, "REFERENCE") == 0)
82 o->opt = CRAM_OPT_REFERENCE, o->val.s = cp;
83 else if (strcmp(arg, "VERSION") == 0)
84 o->opt = CRAM_OPT_VERSION, o->val.s =cp;
85 else if (strcmp(arg, "MULTI_SEQ_PER_SLICE") == 0)
86 o->opt = CRAM_OPT_MULTI_SEQ_PER_SLICE, o->val.i = atoi(cp);
87 else if (strcmp(arg, "NTHREADS") == 0)
88 o->opt = CRAM_OPT_NTHREADS, o->val.i = atoi(cp);
89 else if (strcmp(arg, "REQUIRED_FIELDS") == 0)
90 o->opt = CRAM_OPT_REQUIRED_FIELDS, o->val.i = strtol(cp, NULL, 0);
91 else {
92 fprintf(stderr, "Unknown option '%s'\n", arg);
93 free(o);
94 return -1;
95 }
96
97 o->next = NULL;
98
99 if (*opts) {
100 t = *opts;
101 while (t->next)
102 t = t->next;
103 t->next = o;
104 } else {
105 *opts = o;
106 }
107
108 return 0;
109 }
110
111 int main(int argc, char *argv[])
112 {
113 samFile *in;
114 char *fn_ref = 0;
115 int flag = 0, c, clevel = -1, ignore_sam_err = 0;
116 char moder[8];
117 bam_hdr_t *h;
118 bam1_t *b;
119 htsFile *out;
120 char modew[8];
121 int r = 0, exit_code = 0;
122 hts_opt *in_opts = NULL, *out_opts = NULL, *last = NULL;
123
124 while ((c = getopt(argc, argv, "IbDCSl:t:i:o:")) >= 0) {
125 switch (c) {
126 case 'S': flag |= 1; break;
127 case 'b': flag |= 2; break;
128 case 'D': flag |= 4; break;
129 case 'C': flag |= 8; break;
130 case 'l': clevel = atoi(optarg); flag |= 2; break;
131 case 't': fn_ref = optarg; break;
132 case 'I': ignore_sam_err = 1; break;
133 case 'i': if (add_option(&in_opts, optarg)) return 1; break;
134 case 'o': if (add_option(&out_opts, optarg)) return 1; break;
135 }
136 }
137 if (argc == optind) {
138 fprintf(stderr, "Usage: samview [-bSCSI] [-l level] [-o option=value] <in.bam>|<in.sam>|<in.cram> [region]\n");
139 return 1;
140 }
141 strcpy(moder, "r");
142 if (flag&4) strcat(moder, "c");
143 else if ((flag&1) == 0) strcat(moder, "b");
144
145 in = sam_open(argv[optind], moder);
146 if (in == NULL) {
147 fprintf(stderr, "Error opening \"%s\"\n", argv[optind]);
148 return EXIT_FAILURE;
149 }
150 h = sam_hdr_read(in);
151 h->ignore_sam_err = ignore_sam_err;
152 b = bam_init1();
153
154 strcpy(modew, "w");
155 if (clevel >= 0 && clevel <= 9) sprintf(modew + 1, "%d", clevel);
156 if (flag&8) strcat(modew, "c");
157 else if (flag&2) strcat(modew, "b");
158 out = hts_open("-", modew);
159 if (out == NULL) {
160 fprintf(stderr, "Error opening standard output\n");
161 return EXIT_FAILURE;
162 }
163
164 /* CRAM output */
165 if (flag & 8) {
166 // Parse input header and use for CRAM output
167 out->fp.cram->header = sam_hdr_parse_(h->text, h->l_text);
168
169 // Create CRAM references arrays
170 if (fn_ref)
171 cram_set_option(out->fp.cram, CRAM_OPT_REFERENCE, fn_ref);
172 else
173 // Attempt to fill out a cram->refs[] array from @SQ headers
174 cram_set_option(out->fp.cram, CRAM_OPT_REFERENCE, NULL);
175 }
176
177 // Process any options; currently cram only.
178 for (; in_opts; in_opts = (last=in_opts)->next, free(last)) {
179 hts_set_opt(in, in_opts->opt, in_opts->val);
180 if (in_opts->opt == CRAM_OPT_REFERENCE)
181 hts_set_opt(out, in_opts->opt, in_opts->val);
182 }
183 for (; out_opts; out_opts = (last=out_opts)->next, free(last))
184 hts_set_opt(out, out_opts->opt, out_opts->val);
185
186 sam_hdr_write(out, h);
187 if (optind + 1 < argc && !(flag&1)) { // BAM input and has a region
188 int i;
189 hts_idx_t *idx;
190 if ((idx = bam_index_load(argv[optind])) == 0) {
191 fprintf(stderr, "[E::%s] fail to load the BAM index\n", __func__);
192 return 1;
193 }
194 for (i = optind + 1; i < argc; ++i) {
195 hts_itr_t *iter;
196 if ((iter = bam_itr_querys(idx, h, argv[i])) == 0) {
197 fprintf(stderr, "[E::%s] fail to parse region '%s'\n", __func__, argv[i]);
198 continue;
199 }
200 while ((r = bam_itr_next(in, iter, b)) >= 0) {
201 if (sam_write1(out, h, b) < 0) {
202 fprintf(stderr, "Error writing output.\n");
203 exit_code = 1;
204 break;
205 }
206 }
207 hts_itr_destroy(iter);
208 }
209 hts_idx_destroy(idx);
210 } else while ((r = sam_read1(in, h, b)) >= 0) {
211 if (sam_write1(out, h, b) < 0) {
212 fprintf(stderr, "Error writing output.\n");
213 exit_code = 1;
214 break;
215 }
216 }
217
218 if (r < -1) {
219 fprintf(stderr, "Error parsing input.\n");
220 exit_code = 1;
221 }
222
223 r = sam_close(out);
224 if (r < 0) {
225 fprintf(stderr, "Error closing output.\n");
226 exit_code = 1;
227 }
228
229 bam_destroy1(b);
230 bam_hdr_destroy(h);
231
232 r = sam_close(in);
233 if (r < 0) {
234 fprintf(stderr, "Error closing input.\n");
235 exit_code = 1;
236 }
237
238 return exit_code;
239 }