|
0
|
1 /* test/sam.c -- SAM/BAM/CRAM API test cases.
|
|
|
2
|
|
|
3 Copyright (C) 2014-2015 Genome Research Ltd.
|
|
|
4
|
|
|
5 Author: John Marshall <jm18@sanger.ac.uk>
|
|
|
6
|
|
|
7 Permission is hereby granted, free of charge, to any person obtaining a copy
|
|
|
8 of this software and associated documentation files (the "Software"), to deal
|
|
|
9 in the Software without restriction, including without limitation the rights
|
|
|
10 to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
|
|
|
11 copies of the Software, and to permit persons to whom the Software is
|
|
|
12 furnished to do so, subject to the following conditions:
|
|
|
13
|
|
|
14 The above copyright notice and this permission notice shall be included in
|
|
|
15 all copies or substantial portions of the Software.
|
|
|
16
|
|
|
17 THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
|
|
|
18 IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
|
|
|
19 FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
|
|
|
20 THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
|
|
|
21 LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
|
|
|
22 FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
|
|
|
23 DEALINGS IN THE SOFTWARE. */
|
|
|
24
|
|
|
25 #include <stdarg.h>
|
|
|
26 #include <stdio.h>
|
|
|
27 #include <stdlib.h>
|
|
|
28 #include <string.h>
|
|
|
29 #include <math.h>
|
|
|
30
|
|
|
31 #include "htslib/sam.h"
|
|
|
32 #include "htslib/faidx.h"
|
|
|
33 #include "htslib/kstring.h"
|
|
|
34
|
|
|
35 int status;
|
|
|
36
|
|
|
37 static void fail(const char *fmt, ...)
|
|
|
38 {
|
|
|
39 va_list args;
|
|
|
40
|
|
|
41 fprintf(stderr, "Failed: ");
|
|
|
42 va_start(args, fmt);
|
|
|
43 vfprintf(stderr, fmt, args);
|
|
|
44 va_end(args);
|
|
|
45 fprintf(stderr, "\n");
|
|
|
46
|
|
|
47 status = EXIT_FAILURE;
|
|
|
48 }
|
|
|
49
|
|
|
50 uint8_t *check_bam_aux_get(const bam1_t *aln, const char *tag, char type)
|
|
|
51 {
|
|
|
52 uint8_t *p = bam_aux_get(aln, tag);
|
|
|
53 if (p) {
|
|
|
54 if (*p == type) return p;
|
|
|
55 else fail("%s field of type '%c', expected '%c'\n", tag, *p, type);
|
|
|
56 }
|
|
|
57 else fail("can't find %s field\n", tag);
|
|
|
58
|
|
|
59 return NULL;
|
|
|
60 }
|
|
|
61
|
|
|
62 #define PI 3.141592653589793
|
|
|
63 #define E 2.718281828459045
|
|
|
64 #define HELLO "Hello, world!"
|
|
|
65 #define BEEF "DEADBEEF"
|
|
|
66
|
|
|
67 #define str(x) #x
|
|
|
68 #define xstr(x) str(x)
|
|
|
69
|
|
|
70 static int aux_fields1(void)
|
|
|
71 {
|
|
|
72 static const char sam[] = "data:"
|
|
|
73 "@SQ\tSN:one\tLN:1000\n"
|
|
|
74 "@SQ\tSN:two\tLN:500\n"
|
|
|
75 "r1\t0\tone\t500\t20\t8M\t*\t0\t0\tATGCATGC\tqqqqqqqq\tXA:A:k\tXi:i:37\tXf:f:" xstr(PI) "\tXd:d:" xstr(E) "\tXZ:Z:" HELLO "\tXH:H:" BEEF "\tXB:B:c,-2,0,+2\tZZ:i:1000000\tY1:i:-2147483648\tY2:i:-2147483647\tY3:i:-1\tY4:i:0\tY5:i:1\tY6:i:2147483647\tY7:i:2147483648\tY8:i:4294967295\n";
|
|
|
76
|
|
|
77 // Canonical form of the alignment record above, as output by sam_format1()
|
|
|
78 static const char r1[] = "r1\t0\tone\t500\t20\t8M\t*\t0\t0\tATGCATGC\tqqqqqqqq\tXA:A:k\tXi:i:37\tXf:f:3.14159\tXd:d:2.71828\tXZ:Z:" HELLO "\tXH:H:" BEEF "\tXB:B:c,-2,0,2\tZZ:i:1000000\tY1:i:-2147483648\tY2:i:-2147483647\tY3:i:-1\tY4:i:0\tY5:i:1\tY6:i:2147483647\tY7:i:2147483648\tY8:i:4294967295";
|
|
|
79
|
|
|
80 samFile *in = sam_open(sam, "r");
|
|
|
81 bam_hdr_t *header = sam_hdr_read(in);
|
|
|
82 bam1_t *aln = bam_init1();
|
|
|
83 uint8_t *p;
|
|
|
84 uint32_t n;
|
|
|
85 kstring_t ks = { 0, 0, NULL };
|
|
|
86
|
|
|
87 if (sam_read1(in, header, aln) >= 0) {
|
|
|
88 if ((p = check_bam_aux_get(aln, "XA", 'A')) && bam_aux2A(p) != 'k')
|
|
|
89 fail("XA field is '%c', expected 'k'", bam_aux2A(p));
|
|
|
90
|
|
|
91 if ((p = check_bam_aux_get(aln, "Xi", 'C')) && bam_aux2i(p) != 37)
|
|
|
92 fail("Xi field is %d, expected 37", bam_aux2i(p));
|
|
|
93
|
|
|
94 if ((p = check_bam_aux_get(aln, "Xf", 'f')) && fabs(bam_aux2f(p) - PI) > 1E-6)
|
|
|
95 fail("Xf field is %.12f, expected pi", bam_aux2f(p));
|
|
|
96
|
|
|
97 if ((p = check_bam_aux_get(aln, "Xd", 'd')) && fabs(bam_aux2f(p) - E) > 1E-6)
|
|
|
98 fail("Xf field is %.12f, expected e", bam_aux2f(p));
|
|
|
99
|
|
|
100 if ((p = check_bam_aux_get(aln, "XZ", 'Z')) && strcmp(bam_aux2Z(p), HELLO) != 0)
|
|
|
101 fail("XZ field is \"%s\", expected \"%s\"", bam_aux2Z(p), HELLO);
|
|
|
102
|
|
|
103 if ((p = check_bam_aux_get(aln, "XH", 'H')) && strcmp(bam_aux2Z(p), BEEF) != 0)
|
|
|
104 fail("XH field is \"%s\", expected \"%s\"", bam_aux2Z(p), BEEF);
|
|
|
105
|
|
|
106 // TODO Invent and use bam_aux2B()
|
|
|
107 if ((p = check_bam_aux_get(aln, "XB", 'B')) && ! (memcmp(p, "Bc", 2) == 0 && (memcpy(&n, p+2, 4), n) == 3 && memcmp(p+6, "\xfe\x00\x02", 3) == 0))
|
|
|
108 fail("XB field is %c,..., expected c,-2,0,+2", p[1]);
|
|
|
109
|
|
|
110 if ((p = check_bam_aux_get(aln, "ZZ", 'I')) && bam_aux2i(p) != 1000000)
|
|
|
111 fail("ZZ field is %d, expected 1000000", bam_aux2i(p));
|
|
|
112
|
|
|
113 if ((p = bam_aux_get(aln, "Y1")) && bam_aux2i(p) != -2147483647-1)
|
|
|
114 fail("Y1 field is %d, expected -2^31", bam_aux2i(p));
|
|
|
115
|
|
|
116 if ((p = bam_aux_get(aln, "Y2")) && bam_aux2i(p) != -2147483647)
|
|
|
117 fail("Y2 field is %d, expected -2^31+1", bam_aux2i(p));
|
|
|
118
|
|
|
119 if ((p = bam_aux_get(aln, "Y3")) && bam_aux2i(p) != -1)
|
|
|
120 fail("Y3 field is %d, expected -1", bam_aux2i(p));
|
|
|
121
|
|
|
122 if ((p = bam_aux_get(aln, "Y4")) && bam_aux2i(p) != 0)
|
|
|
123 fail("Y4 field is %d, expected 0", bam_aux2i(p));
|
|
|
124
|
|
|
125 if ((p = bam_aux_get(aln, "Y5")) && bam_aux2i(p) != 1)
|
|
|
126 fail("Y5 field is %d, expected 1", bam_aux2i(p));
|
|
|
127
|
|
|
128 if ((p = bam_aux_get(aln, "Y6")) && bam_aux2i(p) != 2147483647)
|
|
|
129 fail("Y6 field is %d, expected 2^31-1", bam_aux2i(p));
|
|
|
130
|
|
|
131 // TODO Checking these perhaps requires inventing bam_aux2u() or so
|
|
|
132 #if 0
|
|
|
133 if ((p = bam_aux_get(aln, "Y7")) && bam_aux2i(p) != 2147483648)
|
|
|
134 fail("Y7 field is %d, expected 2^31", bam_aux2i(p));
|
|
|
135
|
|
|
136 if ((p = bam_aux_get(aln, "Y8")) && bam_aux2i(p) != 4294967295)
|
|
|
137 fail("Y8 field is %d, expected 2^32-1", bam_aux2i(p));
|
|
|
138 #endif
|
|
|
139
|
|
|
140 if (sam_format1(header, aln, &ks) < 0)
|
|
|
141 fail("can't format record");
|
|
|
142
|
|
|
143 if (strcmp(ks.s, r1) != 0)
|
|
|
144 fail("record formatted incorrectly: \"%s\"", ks.s);
|
|
|
145
|
|
|
146 free(ks.s);
|
|
|
147 }
|
|
|
148 else fail("can't read record");
|
|
|
149
|
|
|
150 bam_destroy1(aln);
|
|
|
151 bam_hdr_destroy(header);
|
|
|
152 sam_close(in);
|
|
|
153
|
|
|
154 return 1;
|
|
|
155 }
|
|
|
156
|
|
|
157 static void iterators1(void)
|
|
|
158 {
|
|
|
159 hts_itr_destroy(sam_itr_queryi(NULL, HTS_IDX_REST, 0, 0));
|
|
|
160 hts_itr_destroy(sam_itr_queryi(NULL, HTS_IDX_NONE, 0, 0));
|
|
|
161 }
|
|
|
162
|
|
|
163 static void faidx1(const char *filename)
|
|
|
164 {
|
|
|
165 int n;
|
|
|
166 faidx_t *fai = fai_load(filename);
|
|
|
167 if (fai == NULL) fail("can't load faidx file");
|
|
|
168
|
|
|
169 n = faidx_fetch_nseq(fai);
|
|
|
170 if (n != 7) fail("faidx_fetch_nseq returned %d, expected 7", n);
|
|
|
171
|
|
|
172 n = faidx_nseq(fai);
|
|
|
173 if (n != 7) fail("faidx_nseq returned %d, expected 7", n);
|
|
|
174
|
|
|
175 fai_destroy(fai);
|
|
|
176 }
|
|
|
177
|
|
|
178 int main(int argc, char **argv)
|
|
|
179 {
|
|
|
180 status = EXIT_SUCCESS;
|
|
|
181
|
|
|
182 aux_fields1();
|
|
|
183 iterators1();
|
|
|
184 if (argc >= 2) faidx1(argv[1]);
|
|
|
185
|
|
|
186 return status;
|
|
|
187 }
|
