Mercurial > repos > youngkim > ezbamqc

comparison ezBAMQC/src/htslib/test/sam.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/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 }