|
0
|
1 /* The MIT License
|
|
|
2
|
|
|
3 Copyright (c) 2011 Heng Li <lh3@live.co.uk>
|
|
|
4
|
|
|
5 Permission is hereby granted, free of charge, to any person obtaining
|
|
|
6 a copy of this software and associated documentation files (the
|
|
|
7 "Software"), to deal in the Software without restriction, including
|
|
|
8 without limitation the rights to use, copy, modify, merge, publish,
|
|
|
9 distribute, sublicense, and/or sell copies of the Software, and to
|
|
|
10 permit persons to whom the Software is furnished to do so, subject to
|
|
|
11 the following conditions:
|
|
|
12
|
|
|
13 The above copyright notice and this permission notice shall be
|
|
|
14 included in all copies or substantial portions of the Software.
|
|
|
15
|
|
|
16 THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
|
|
|
17 EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
|
|
|
18 MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
|
|
|
19 NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS
|
|
|
20 BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
|
|
|
21 ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
|
|
|
22 CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
|
|
|
23 SOFTWARE.
|
|
|
24 */
|
|
|
25
|
|
|
26 #include <stdio.h>
|
|
|
27 #include <stdlib.h>
|
|
|
28 #include <string.h>
|
|
|
29 #include <unistd.h>
|
|
|
30 #include <zlib.h>
|
|
|
31 #include "kstring.h"
|
|
|
32 #include "kseq.h"
|
|
|
33 KSTREAM_INIT(gzFile, gzread, 16384)
|
|
|
34
|
|
|
35 #define N_TMPSTR 5
|
|
|
36 #define LINE_LEN 60
|
|
|
37
|
|
|
38 // append a CIGAR operation plus length
|
|
|
39 #define write_cigar(_c, _n, _m, _v) do { \
|
|
|
40 if (_n == _m) { \
|
|
|
41 _m = _m? _m<<1 : 4; \
|
|
|
42 _c = realloc(_c, _m * sizeof(unsigned)); \
|
|
|
43 } \
|
|
|
44 _c[_n++] = (_v); \
|
|
|
45 } while (0)
|
|
|
46
|
|
|
47 // a fatal error
|
|
|
48 static void fatal(const char *msg)
|
|
|
49 {
|
|
|
50 fprintf(stderr, "E %s\n", msg);
|
|
|
51 exit(1);
|
|
|
52 }
|
|
|
53 // remove pads
|
|
|
54 static void remove_pads(const kstring_t *src, kstring_t *dst)
|
|
|
55 {
|
|
|
56 int i, j;
|
|
|
57 dst->l = 0;
|
|
|
58 kputsn(src->s, src->l, dst);
|
|
|
59 for (i = j = 0; i < dst->l; ++i)
|
|
|
60 if (dst->s[i] != '*') dst->s[j++] = dst->s[i];
|
|
|
61 dst->s[j] = 0;
|
|
|
62 dst->l = j;
|
|
|
63 }
|
|
|
64
|
|
|
65 int main(int argc, char *argv[])
|
|
|
66 {
|
|
|
67 gzFile fp;
|
|
|
68 kstream_t *ks;
|
|
|
69 kstring_t s, t[N_TMPSTR];
|
|
|
70 int dret, i, k, af_n, af_max, af_i, c, is_padded = 0, write_cns = 0, *p2u = 0;
|
|
|
71 long m_cigar = 0, n_cigar = 0;
|
|
|
72 unsigned *af, *cigar = 0;
|
|
|
73
|
|
|
74 while ((c = getopt(argc, argv, "pc")) >= 0) {
|
|
|
75 switch (c) {
|
|
|
76 case 'p': is_padded = 1; break;
|
|
|
77 case 'c': write_cns = 1; break;
|
|
|
78 }
|
|
|
79 }
|
|
|
80 if (argc == optind) {
|
|
|
81 fprintf(stderr, "\nUsage: ace2sam [-pc] <in.ace>\n\n");
|
|
|
82 fprintf(stderr, "Options: -p output padded SAM\n");
|
|
|
83 fprintf(stderr, " -c write the contig sequence in SAM\n\n");
|
|
|
84 fprintf(stderr, "Notes: 1. Fields must appear in the following order: (CO->[BQ]->(AF)->(RD->QA))\n");
|
|
|
85 fprintf(stderr, " 2. The order of reads in AF and in RD must be identical\n");
|
|
|
86 fprintf(stderr, " 3. Except in BQ, words and numbers must be separated by a single SPACE or TAB\n");
|
|
|
87 fprintf(stderr, " 4. This program writes the headerless SAM to stdout and header to stderr\n\n");
|
|
|
88 return 1;
|
|
|
89 }
|
|
|
90
|
|
|
91 s.l = s.m = 0; s.s = 0;
|
|
|
92 af_n = af_max = af_i = 0; af = 0;
|
|
|
93 for (i = 0; i < N_TMPSTR; ++i) t[i].l = t[i].m = 0, t[i].s = 0;
|
|
|
94 fp = strcmp(argv[1], "-")? gzopen(argv[optind], "r") : gzdopen(fileno(stdin), "r");
|
|
|
95 ks = ks_init(fp);
|
|
|
96 while (ks_getuntil(ks, 0, &s, &dret) >= 0) {
|
|
|
97 if (strcmp(s.s, "CO") == 0) { // contig sequence
|
|
|
98 kstring_t *cns;
|
|
|
99 t[0].l = t[1].l = t[2].l = t[3].l = t[4].l = 0; // 0: name; 1: padded ctg; 2: unpadded ctg/padded read; 3: unpadded read; 4: SAM line
|
|
|
100 af_n = af_i = 0; // reset the af array
|
|
|
101 ks_getuntil(ks, 0, &s, &dret); kputs(s.s, &t[0]); // contig name
|
|
|
102 ks_getuntil(ks, '\n', &s, &dret); // read the whole line
|
|
|
103 while (ks_getuntil(ks, '\n', &s, &dret) >= 0 && s.l > 0) kputsn(s.s, s.l, &t[1]); // read the padded consensus sequence
|
|
|
104 remove_pads(&t[1], &t[2]); // construct the unpadded sequence
|
|
|
105 // compute the array for mapping padded positions to unpadded positions
|
|
|
106 p2u = realloc(p2u, t[1].m * sizeof(int));
|
|
|
107 for (i = k = 0; i < t[1].l; ++i) {
|
|
|
108 p2u[i] = k;
|
|
|
109 if (t[1].s[i] != '*') ++k;
|
|
|
110 }
|
|
|
111 // write out the SAM header and contig sequences
|
|
|
112 fprintf(stderr, "H @SQ\tSN:%s\tLN:%ld\n", t[0].s, t[is_padded?1:2].l); // The SAM header line
|
|
|
113 cns = &t[is_padded?1:2];
|
|
|
114 fprintf(stderr, "S >%s\n", t[0].s);
|
|
|
115 for (i = 0; i < cns->l; i += LINE_LEN) {
|
|
|
116 fputs("S ", stderr);
|
|
|
117 for (k = 0; k < LINE_LEN && i + k < cns->l; ++k)
|
|
|
118 fputc(cns->s[i + k], stderr);
|
|
|
119 fputc('\n', stderr);
|
|
|
120 }
|
|
|
121
|
|
|
122 #define __padded2cigar(sp) do { \
|
|
|
123 int i, l_M = 0, l_D = 0; \
|
|
|
124 for (i = 0; i < sp.l; ++i) { \
|
|
|
125 if (sp.s[i] == '*') { \
|
|
|
126 if (l_M) write_cigar(cigar, n_cigar, m_cigar, l_M<<4); \
|
|
|
127 ++l_D; l_M = 0; \
|
|
|
128 } else { \
|
|
|
129 if (l_D) write_cigar(cigar, n_cigar, m_cigar, l_D<<4 | 2); \
|
|
|
130 ++l_M; l_D = 0; \
|
|
|
131 } \
|
|
|
132 } \
|
|
|
133 if (l_M) write_cigar(cigar, n_cigar, m_cigar, l_M<<4); \
|
|
|
134 else write_cigar(cigar, n_cigar, m_cigar, l_D<<4 | 2); \
|
|
|
135 } while (0)
|
|
|
136
|
|
|
137 if (write_cns) { // write the consensus SAM line (dummy read)
|
|
|
138 n_cigar = 0;
|
|
|
139 if (is_padded) __padded2cigar(t[1]);
|
|
|
140 else write_cigar(cigar, n_cigar, m_cigar, t[2].l<<4);
|
|
|
141 kputsn(t[0].s, t[0].l, &t[4]); kputs("\t516\t", &t[4]); kputsn(t[0].s, t[0].l, &t[4]); kputs("\t1\t60\t", &t[4]);
|
|
|
142 for (i = 0; i < n_cigar; ++i) {
|
|
|
143 kputw(cigar[i]>>4, &t[4]); kputc("MIDNSHP=X"[cigar[i]&0xf], &t[4]);
|
|
|
144 }
|
|
|
145 kputs("\t*\t0\t0\t", &t[4]); kputsn(t[2].s, t[2].l, &t[4]); kputs("\t*", &t[4]);
|
|
|
146 }
|
|
|
147 } else if (strcmp(s.s, "BQ") == 0) { // contig quality
|
|
|
148 if (t[0].l == 0) fatal("come to 'BQ' before reading 'CO'");
|
|
|
149 if (dret != '\n') ks_getuntil(ks, '\n', &s, &dret); // read the entire "BQ" line
|
|
|
150 if (write_cns) t[4].s[--t[4].l] = 0; // remove the trailing "*"
|
|
|
151 for (i = 0; i < t[2].l; ++i) { // read the consensus quality
|
|
|
152 int q;
|
|
|
153 if (ks_getuntil(ks, 0, &s, &dret) < 0) fprintf(stderr, "E truncated contig quality\n");
|
|
|
154 if (s.l) {
|
|
|
155 q = atoi(s.s) + 33;
|
|
|
156 if (q > 126) q = 126;
|
|
|
157 if (write_cns) kputc(q, &t[4]);
|
|
|
158 } else --i;
|
|
|
159 }
|
|
|
160 if (dret != '\n') ks_getuntil(ks, '\n', &s, &dret);
|
|
|
161 ks_getuntil(ks, '\n', &s, &dret); // skip the empty line
|
|
|
162 if (write_cns) puts(t[4].s); t[4].l = 0;
|
|
|
163 } else if (strcmp(s.s, "AF") == 0) { // padded read position
|
|
|
164 int reversed, neg, pos;
|
|
|
165 if (t[0].l == 0) fatal("come to 'AF' before reading 'CO'");
|
|
|
166 if (write_cns) {
|
|
|
167 if (t[4].l) puts(t[4].s);
|
|
|
168 t[4].l = 0;
|
|
|
169 }
|
|
|
170 ks_getuntil(ks, 0, &s, &dret); // read name
|
|
|
171 ks_getuntil(ks, 0, &s, &dret); reversed = s.s[0] == 'C'? 1 : 0; // strand
|
|
|
172 ks_getuntil(ks, 0, &s, &dret); pos = atoi(s.s); neg = pos < 0? 1 : 0; pos = pos < 0? -pos : pos; // position
|
|
|
173 if (af_n == af_max) { // double the af array
|
|
|
174 af_max = af_max? af_max<<1 : 4;
|
|
|
175 af = realloc(af, af_max * sizeof(unsigned));
|
|
|
176 }
|
|
|
177 af[af_n++] = pos << 2 | neg << 1 | reversed; // keep the placement information
|
|
|
178 } else if (strcmp(s.s, "RD") == 0) { // read sequence
|
|
|
179 if (af_i >= af_n) fatal("more 'RD' records than 'AF'");
|
|
|
180 t[2].l = t[3].l = t[4].l = 0;
|
|
|
181 ks_getuntil(ks, 0, &t[4], &dret); // QNAME
|
|
|
182 if (dret != '\n') ks_getuntil(ks, '\n', &s, &dret); // read the entire RD line
|
|
|
183 while (ks_getuntil(ks, '\n', &s, &dret) >= 0 && s.l > 0) kputs(s.s, &t[2]); // read the read sequence
|
|
|
184 } else if (strcmp(s.s, "QA") == 0) { // clipping
|
|
|
185 if (af_i >= af_n) fatal("more 'QA' records than 'AF'");
|
|
|
186 int beg, end, pos, op;
|
|
|
187 ks_getuntil(ks, 0, &s, &dret); ks_getuntil(ks, 0, &s, &dret); // skip quality clipping
|
|
|
188 ks_getuntil(ks, 0, &s, &dret); beg = atoi(s.s) - 1; // align clipping start
|
|
|
189 ks_getuntil(ks, 0, &s, &dret); end = atoi(s.s); // clipping end
|
|
|
190 if (dret != '\n') ks_getuntil(ks, '\n', &s, &dret);
|
|
|
191 // compute 1-based POS
|
|
|
192 pos = af[af_i]>>2; // retrieve the position information
|
|
|
193 if (af[af_i]>>1&1) pos = -pos;
|
|
|
194 pos += beg; // now pos is the true padded position
|
|
|
195 // generate CIGAR
|
|
|
196 remove_pads(&t[2], &t[3]); // backup the unpadded read sequence
|
|
|
197 n_cigar = 0;
|
|
|
198 if (beg) write_cigar(cigar, n_cigar, m_cigar, beg<<4|4);
|
|
|
199 if (is_padded) {
|
|
|
200 __padded2cigar(t[2]);
|
|
|
201 if (beg && n_cigar > 1) cigar[1] -= beg<<4; // fix the left-hand CIGAR
|
|
|
202 if (end < t[2].l && n_cigar) cigar[n_cigar-1] -= (t[2].l - end)<<4; // fix the right-hand CIGAR
|
|
|
203 } else {
|
|
|
204 // generate flattened CIGAR string
|
|
|
205 for (i = beg, k = pos - 1; i < end; ++i, ++k)
|
|
|
206 t[2].s[i] = t[2].s[i] != '*'? (t[1].s[k] != '*'? 0 : 1) : (t[1].s[k] != '*'? 2 : 6);
|
|
|
207 // generate the proper CIGAR
|
|
|
208 for (i = beg + 1, k = 1, op = t[2].s[beg]; i < end; ++i) {
|
|
|
209 if (op != t[2].s[i]) {
|
|
|
210 write_cigar(cigar, n_cigar, m_cigar, k<<4|op);
|
|
|
211 op = t[2].s[i]; k = 1;
|
|
|
212 } else ++k;
|
|
|
213 }
|
|
|
214 write_cigar(cigar, n_cigar, m_cigar, k<<4|op);
|
|
|
215 // remove unnecessary "P" and possibly merge adjacent operations
|
|
|
216 for (i = 2; i < n_cigar; ++i) {
|
|
|
217 if ((cigar[i]&0xf) != 1 && (cigar[i-1]&0xf) == 6 && (cigar[i-2]&0xf) != 1) {
|
|
|
218 cigar[i-1] = 0;
|
|
|
219 if ((cigar[i]&0xf) == (cigar[i-2]&0xf)) // merge operations
|
|
|
220 cigar[i] += cigar[i-2], cigar[i-2] = 0;
|
|
|
221 }
|
|
|
222 }
|
|
|
223 for (i = k = 0; i < n_cigar; ++i) // squeeze out dumb operations
|
|
|
224 if (cigar[i]) cigar[k++] = cigar[i];
|
|
|
225 n_cigar = k;
|
|
|
226 }
|
|
|
227 if (end < t[2].l) write_cigar(cigar, n_cigar, m_cigar, (t[2].l - end)<<4|4);
|
|
|
228 // write the SAM line for the read
|
|
|
229 kputc('\t', &t[4]); // QNAME has already been written
|
|
|
230 kputw((af[af_i]&1)? 16 : 0, &t[4]); kputc('\t', &t[4]); // FLAG
|
|
|
231 kputsn(t[0].s, t[0].l, &t[4]); kputc('\t', &t[4]); // RNAME
|
|
|
232 kputw(is_padded? pos : p2u[pos-1]+1, &t[4]); // POS
|
|
|
233 kputs("\t60\t", &t[4]); // MAPQ
|
|
|
234 for (i = 0; i < n_cigar; ++i) { // CIGAR
|
|
|
235 kputw(cigar[i]>>4, &t[4]); kputc("MIDNSHP=X"[cigar[i]&0xf], &t[4]);
|
|
|
236 }
|
|
|
237 kputs("\t*\t0\t0\t", &t[4]); // empty MRNM, MPOS and TLEN
|
|
|
238 kputsn(t[3].s, t[3].l, &t[4]); // unpadded SEQ
|
|
|
239 kputs("\t*", &t[4]); // QUAL
|
|
|
240 puts(t[4].s); // print to stdout
|
|
|
241 ++af_i;
|
|
|
242 } else if (dret != '\n') ks_getuntil(ks, '\n', &s, &dret);
|
|
|
243 }
|
|
|
244 ks_destroy(ks);
|
|
|
245 gzclose(fp);
|
|
|
246 free(af); free(s.s); free(cigar); free(p2u);
|
|
|
247 for (i = 0; i < N_TMPSTR; ++i) free(t[i].s);
|
|
|
248 return 0;
|
|
|
249 }
|
