1 |
#ifndef _G_BAM_H |
2 |
#define _G_BAM_H |
3 |
#include "GBase.h" |
4 |
#include "GList.hh" |
5 |
#include "bam.h" |
6 |
#include "sam.h" |
7 |
|
8 |
|
9 |
class GBamReader; |
10 |
class GBamWriter; |
11 |
|
12 |
class GBamRecord: public GSeg { |
13 |
friend class GBamReader; |
14 |
friend class GBamWriter; |
15 |
bam1_t* b; |
16 |
// b->data has the following strings concatenated: |
17 |
// qname (including the terminal \0) |
18 |
// +cigar (each event encoded on 32 bits) |
19 |
// +seq (4bit-encoded) |
20 |
// +qual |
21 |
// +aux |
22 |
//bool novel; //created from scratch, not from a given bam1_t* record |
23 |
bam_header_t* bam_header; |
24 |
|
25 |
public: |
26 |
GVec<GSeg> exons; |
27 |
//created from a reader: |
28 |
GBamRecord(bam1_t* from_b=NULL, bam_header_t* b_header=NULL):exons(1) { |
29 |
bam_header=NULL; |
30 |
if (from_b==NULL) { |
31 |
b=bam_init1(); |
32 |
//novel=true; |
33 |
} |
34 |
else { |
35 |
b=from_b; //it'll take over from_b and FREE it when destroyed |
36 |
//novel=false; |
37 |
} |
38 |
bam_header=b_header; |
39 |
setupCoordinates();//sets start, end coordinate and populates exons array |
40 |
} |
41 |
|
42 |
const GBamRecord& operator=(GBamRecord& r) { |
43 |
//make a new copy of the bam1_t record etc. |
44 |
clear(); |
45 |
b=bam_dup1(r.b); |
46 |
//novel=true; //will also free b when destroyed |
47 |
start=r.start; |
48 |
end=r.end; |
49 |
exons = r.exons; |
50 |
return *this; |
51 |
} |
52 |
|
53 |
void setupCoordinates(); |
54 |
void clear() { |
55 |
//if (novel) { |
56 |
bam_destroy1(b); |
57 |
// novel=false; |
58 |
// } |
59 |
b=NULL; |
60 |
exons.Clear(); |
61 |
bam_header=NULL; |
62 |
} |
63 |
|
64 |
~GBamRecord() { |
65 |
clear(); |
66 |
} |
67 |
|
68 |
void parse_error(const char* s) { |
69 |
GError("BAM parsing error: %s\n", s); |
70 |
} |
71 |
|
72 |
bam1_t* get_b() { return b; } |
73 |
|
74 |
void set_mdata(int32_t mtid, int32_t m0pos, //0-based coordinate, -1 if not available |
75 |
int32_t isize=0) { //mate info for current record |
76 |
b->core.mtid=mtid; |
77 |
b->core.mpos=m0pos; // should be -1 if '*' |
78 |
b->core.isize=isize; //should be 0 if not available |
79 |
} |
80 |
|
81 |
void set_flags(uint16_t flags) { |
82 |
b->core.flag=flags; |
83 |
} |
84 |
|
85 |
//creates a new record from 1-based alignment coordinate |
86 |
//quals should be given as Phred33 |
87 |
//Warning: pos and mate_pos must be given 1-based! |
88 |
GBamRecord(const char* qname, int32_t gseq_tid, |
89 |
int pos, bool reverse, const char* qseq, const char* cigar=NULL, const char* quals=NULL); |
90 |
GBamRecord(const char* qname, int32_t flags, int32_t g_tid, |
91 |
int pos, int map_qual, const char* cigar, int32_t mg_tid, int mate_pos, |
92 |
int insert_size, const char* qseq, const char* quals=NULL, |
93 |
GVec<char*>* aux_strings=NULL); |
94 |
//const std::vector<std::string>* aux_strings=NULL); |
95 |
void set_cigar(const char* cigar); //converts and adds CIGAR string given in plain SAM text format |
96 |
void add_sequence(const char* qseq, int slen=-1); //adds the DNA sequence given in plain text format |
97 |
void add_quals(const char* quals); //quality values string in Phred33 format |
98 |
void add_aux(const char* str); //adds one aux field in plain SAM text format (e.g. "NM:i:1") |
99 |
void add_aux(const char tag[2], char atype, int len, uint8_t *data) { |
100 |
int ori_len = b->data_len; |
101 |
b->data_len += 3 + len; |
102 |
b->l_aux += 3 + len; |
103 |
if (b->m_data < b->data_len) { |
104 |
b->m_data = b->data_len; |
105 |
kroundup32(b->m_data); |
106 |
b->data = (uint8_t*)realloc(b->data, b->m_data); |
107 |
} |
108 |
b->data[ori_len] = tag[0]; b->data[ori_len + 1] = tag[1]; |
109 |
b->data[ori_len + 2] = atype; |
110 |
memcpy(b->data + ori_len + 3, data, len); |
111 |
} |
112 |
void add_tag(const char tag[2], char atype, int len, uint8_t *data) { |
113 |
//same with add_aux() |
114 |
add_aux(tag,atype,len,data); |
115 |
} |
116 |
//--query methods: |
117 |
uint32_t flags() { return b->core.flag; } |
118 |
bool isUnmapped() { return ((b->core.flag & BAM_FUNMAP) != 0); } |
119 |
bool isMapped() { return ((b->core.flag & BAM_FUNMAP) == 0); } |
120 |
bool isPaired() { return ((b->core.flag & BAM_FPAIRED) != 0); } |
121 |
const char* name() { return bam1_qname(b); } |
122 |
int pairOrder() { |
123 |
//which read in the pair: 0 = unpaired, 1=first read, 2=second read |
124 |
int r=0; |
125 |
if ((b->core.flag & BAM_FREAD1) != 0) r=1; |
126 |
else if ((b->core.flag & BAM_FREAD2) != 0) r=2; |
127 |
return r; |
128 |
} |
129 |
bool revStrand() { |
130 |
//this is the raw alignment strand, NOT the transcription (XS) strand |
131 |
return ((b->core.flag & BAM_FREVERSE) != 0); |
132 |
} |
133 |
const char* refName() { |
134 |
return (bam_header!=NULL) ? |
135 |
((b->core.tid<0) ? "*" : bam_header->target_name[b->core.tid]) : NULL; |
136 |
} |
137 |
int32_t refId() { return b->core.tid; } |
138 |
int32_t mate_refId() { return b->core.mtid; } |
139 |
const char* mate_refName() { |
140 |
return (bam_header!=NULL) ? |
141 |
((b->core.mtid<0) ? "*" : bam_header->target_name[b->core.mtid]) : NULL; |
142 |
} |
143 |
int32_t insertSize() { return b->core.isize; } |
144 |
int32_t mate_start() { return b->core.mpos<0? 0 : b->core.mpos+1; } |
145 |
|
146 |
int find_tag(const char tag[2], uint8_t* & s, char& tag_type); |
147 |
//position s at the beginning of tag data, tag_type is set to the found tag type |
148 |
//returns length of tag data, or 0 if tag not found |
149 |
|
150 |
char* tag_str(const char tag[2]); //return tag value for tag type 'Z' |
151 |
int tag_int(const char tag[2]); //return numeric value of tag (for numeric types) |
152 |
char tag_char(const char tag[2]); //return char value of tag (for type 'A') |
153 |
char spliceStrand(); // '+', '-' from the XS tag, or '.' if no XS tag |
154 |
|
155 |
char* sequence(); //user should free after use |
156 |
char* qualities();//user should free after use |
157 |
char* cigar(); //returns text version of the CIGAR string; user must free |
158 |
}; |
159 |
|
160 |
|
161 |
class GBamReader { |
162 |
samfile_t* bam_file; |
163 |
char* fname; |
164 |
public: |
165 |
void bopen(const char* filename) { |
166 |
if (strcmp(filename, "-")==0) { |
167 |
//if stdin was given, we HAVE to assume it's text SAM format, sorry |
168 |
bam_file=samopen(filename, "r", 0); |
169 |
} |
170 |
else { |
171 |
//BAM files have the zlib signature bytes at the beginning: 1F 8B 08 |
172 |
//if that's not present, we assume sam file |
173 |
FILE* f=fopen(filename, "rb"); |
174 |
if (f==NULL) GError("Error opening file %s!\n", filename); |
175 |
byte fsig[3]; |
176 |
size_t rd=fread(fsig, 1, 3, f); |
177 |
if (rd<3) GError("Error reading from file %s!\n",filename); |
178 |
if ((fsig[0]==0x1F && fsig[1]==0x8B && fsig[2]==0x08) || |
179 |
(fsig[0]=='B' && fsig[1]=='A' && fsig[2]=='M')) { |
180 |
bam_file=samopen(filename, "rb", 0); //BAM or uncompressed BAM |
181 |
} |
182 |
else { //assume text SAM file |
183 |
bam_file=samopen(filename, "r", 0); |
184 |
} |
185 |
} |
186 |
if (bam_file==NULL) |
187 |
GError("Error: could not open SAM file %s!\n",filename); |
188 |
fname=Gstrdup(filename); |
189 |
} |
190 |
GBamReader(const char* fn) { |
191 |
bam_file=NULL; |
192 |
fname=NULL; |
193 |
bopen(fn); |
194 |
} |
195 |
|
196 |
bam_header_t* header() { |
197 |
return bam_file? bam_file->header : NULL; |
198 |
} |
199 |
void bclose() { |
200 |
if (bam_file) { |
201 |
samclose(bam_file); |
202 |
bam_file=NULL; |
203 |
} |
204 |
} |
205 |
~GBamReader() { |
206 |
bclose(); |
207 |
GFREE(fname); |
208 |
} |
209 |
void rewind() { |
210 |
if (fname==NULL) { |
211 |
GMessage("Warning: GBamReader::rewind() called without a file name.\n"); |
212 |
return; |
213 |
} |
214 |
bclose(); |
215 |
char* ifname=fname; |
216 |
bopen(ifname); |
217 |
GFREE(ifname); |
218 |
} |
219 |
|
220 |
GBamRecord* next() { |
221 |
if (bam_file==NULL) |
222 |
GError("Warning: GBamReader::next() called with no open file.\n"); |
223 |
bam1_t* b = bam_init1(); |
224 |
if (samread(bam_file, b) >= 0) { |
225 |
GBamRecord* bamrec=new GBamRecord(b, bam_file->header); |
226 |
return bamrec; |
227 |
} |
228 |
else { |
229 |
bam_destroy1(b); |
230 |
return NULL; |
231 |
} |
232 |
} |
233 |
}; |
234 |
|
235 |
|
236 |
class GBamWriter { |
237 |
samfile_t* bam_file; |
238 |
bam_header_t* bam_header; |
239 |
public: |
240 |
void create(const char* fname, bool uncompressed=false) { |
241 |
if (bam_header==NULL) |
242 |
GError("Error: no bam_header for GBamWriter::create()!\n"); |
243 |
if (uncompressed) { |
244 |
bam_file=samopen(fname, "wbu", bam_header); |
245 |
} |
246 |
else { |
247 |
bam_file=samopen(fname, "wb", bam_header); |
248 |
} |
249 |
if (bam_file==NULL) |
250 |
GError("Error: could not create BAM file %s!\n",fname); |
251 |
} |
252 |
void create(const char* fname, bam_header_t* bh, bool uncompressed=false) { |
253 |
bam_header=bh; |
254 |
create(fname,uncompressed); |
255 |
} |
256 |
|
257 |
GBamWriter(const char* fname, bam_header_t* bh, bool uncompressed=false) { |
258 |
create(fname, bh, uncompressed); |
259 |
} |
260 |
|
261 |
GBamWriter(const char* fname, const char* samfname, bool uncompressed=false) { |
262 |
tamFile samf_in=sam_open(samfname); |
263 |
if (samf_in==NULL) |
264 |
GError("Error: could not open SAM file %s\n", samfname); |
265 |
bam_header=sam_header_read(samf_in); |
266 |
if (bam_header==NULL) |
267 |
GError("Error: could not read SAM header from %s!\n",samfname); |
268 |
sam_close(samf_in); |
269 |
create(fname, uncompressed); |
270 |
} |
271 |
|
272 |
~GBamWriter() { |
273 |
samclose(bam_file); |
274 |
bam_header_destroy(bam_header); |
275 |
} |
276 |
bam_header_t* get_header() { return bam_header; } |
277 |
int32_t get_tid(const char *seq_name) { |
278 |
if (bam_header==NULL) |
279 |
GError("Error: missing SAM header (get_tid())\n"); |
280 |
return bam_get_tid(bam_header, seq_name); |
281 |
} |
282 |
|
283 |
//just a convenience function for creating a new record, but it's NOT written |
284 |
//given pos must be 1-based (so it'll be stored as pos-1 because BAM is 0-based) |
285 |
GBamRecord* new_record(const char* qname, const char* gseqname, |
286 |
int pos, bool reverse, const char* qseq, const char* cigar=NULL, const char* qual=NULL) { |
287 |
int32_t gseq_tid=get_tid(gseqname); |
288 |
if (gseq_tid < 0 && strcmp(gseqname, "*")) { |
289 |
if (bam_header->n_targets == 0) { |
290 |
GError("Error: missing/invalid SAM header\n"); |
291 |
} else |
292 |
GMessage("Warning: reference '%s' not found in header, will consider it '*'.\n", |
293 |
gseqname); |
294 |
} |
295 |
|
296 |
return (new GBamRecord(qname, gseq_tid, pos, reverse, qseq, cigar, qual)); |
297 |
} |
298 |
|
299 |
GBamRecord* new_record(const char* qname, int32_t flags, const char* gseqname, |
300 |
int pos, int map_qual, const char* cigar, const char* mgseqname, int mate_pos, |
301 |
int insert_size, const char* qseq, const char* quals=NULL, |
302 |
GVec<char*>* aux_strings=NULL) { |
303 |
int32_t gseq_tid=get_tid(gseqname); |
304 |
if (gseq_tid < 0 && strcmp(gseqname, "*")) { |
305 |
if (bam_header->n_targets == 0) { |
306 |
GError("Error: missing/invalid SAM header\n"); |
307 |
} else |
308 |
GMessage("Warning: reference '%s' not found in header, will consider it '*'.\n", |
309 |
gseqname); |
310 |
} |
311 |
int32_t mgseq_tid=-1; |
312 |
if (mgseqname!=NULL) { |
313 |
if (strcmp(mgseqname, "=")==0) { |
314 |
mgseq_tid=gseq_tid; |
315 |
} |
316 |
else { |
317 |
mgseq_tid=get_tid(mgseqname); |
318 |
if (mgseq_tid < 0 && strcmp(mgseqname, "*")) { |
319 |
GMessage("Warning: reference '%s' not found in header, will consider it '*'.\n", |
320 |
mgseqname); |
321 |
} |
322 |
} |
323 |
} |
324 |
return (new GBamRecord(qname, flags, gseq_tid, pos, map_qual, cigar, |
325 |
mgseq_tid, mate_pos, insert_size, qseq, quals, aux_strings)); |
326 |
} |
327 |
|
328 |
void write(GBamRecord* brec) { |
329 |
if (brec!=NULL) |
330 |
samwrite(this->bam_file,brec->get_b()); |
331 |
} |
332 |
void write(bam1_t* b) { |
333 |
samwrite(this->bam_file, b); |
334 |
} |
335 |
}; |
336 |
|
337 |
#endif |