gclib/gclib/GBam.h

#ifndef _G_BAM_H
#define _G_BAM_H
#include "GBase.h"
#include "GList.hh"
#include "bam.h"
#include "sam.h"


class GBamReader;
class GBamWriter;

class GBamRecord: public GSeg {
   friend class GBamReader;
   friend class GBamWriter;
   bam1_t* b;
   // b->data has the following strings concatenated:
   //  qname (including the terminal \0)
   //  +cigar (each event encoded on 32 bits)
   //   +seq  (4bit-encoded)
   //    +qual
   //     +aux
   //bool novel; //created from scratch, not from a given bam1_t* record
   bam_header_t* bam_header;

 public:
   GVec<GSeg> exons;
   //created from a reader:
   GBamRecord(bam1_t* from_b=NULL, bam_header_t* b_header=NULL):exons(1) {
      bam_header=NULL;
      if (from_b==NULL) {
           b=bam_init1();
           //novel=true;
           }
        else {
           b=from_b; //it'll take over from_b and FREE it when destroyed
           //novel=false;
           }
      bam_header=b_header;
      setupCoordinates();//sets start, end coordinate and populates exons array
      }

   const GBamRecord& operator=(GBamRecord& r) {
      //make a new copy of the bam1_t record etc.
      clear();
      b=bam_dup1(r.b);
      //novel=true; //will also free b when destroyed
      start=r.start;
      end=r.end;
      exons = r.exons;
      return *this;
      }

     void setupCoordinates();
     void clear() {
        //if (novel) {
        bam_destroy1(b);
        //   novel=false;
        //   }
        b=NULL;
        exons.Clear();
        bam_header=NULL;
        }

    ~GBamRecord() {
       clear();
       }

    void parse_error(const char* s) {
      GError("BAM parsing error: %s\n", s);
      }

    bam1_t* get_b() { return b; }

    void set_mdata(int32_t mtid, int32_t m0pos, //0-based coordinate, -1 if not available
                     int32_t isize=0) { //mate info for current record
      b->core.mtid=mtid;
      b->core.mpos=m0pos; // should be -1 if '*'
      b->core.isize=isize; //should be 0 if not available
      }

    void set_flags(uint16_t flags) {
      b->core.flag=flags;
      }

    //creates a new record from 1-based alignment coordinate
    //quals should be given as Phred33
    //Warning: pos and mate_pos must be given 1-based!
    GBamRecord(const char* qname, int32_t gseq_tid,
                    int pos, bool reverse, const char* qseq, const char* cigar=NULL, const char* quals=NULL);
    GBamRecord(const char* qname, int32_t flags, int32_t g_tid,
             int pos, int map_qual, const char* cigar, int32_t mg_tid, int mate_pos,
             int insert_size, const char* qseq, const char* quals=NULL,
             GVec<char*>* aux_strings=NULL);
             //const std::vector<std::string>* aux_strings=NULL);
    void set_cigar(const char* cigar); //converts and adds CIGAR string given in plain SAM text format
    void add_sequence(const char* qseq, int slen=-1); //adds the DNA sequence given in plain text format
    void add_quals(const char* quals); //quality values string in Phred33 format
    void add_aux(const char* str); //adds one aux field in plain SAM text format (e.g. "NM:i:1")
    void add_aux(const char tag[2], char atype, int len, uint8_t *data) {
      int ori_len = b->data_len;
      b->data_len += 3 + len;
      b->l_aux += 3 + len;
      if (b->m_data < b->data_len) {
        b->m_data = b->data_len;
        kroundup32(b->m_data);
        b->data = (uint8_t*)realloc(b->data, b->m_data);
      }
      b->data[ori_len] = tag[0]; b->data[ori_len + 1] = tag[1];
      b->data[ori_len + 2] = atype;
      memcpy(b->data + ori_len + 3, data, len);
      }
    void add_tag(const char tag[2], char atype, int len, uint8_t *data) {
      //same with add_aux()
      add_aux(tag,atype,len,data);
      }
 //--query methods:
 uint32_t flags() { return b->core.flag; }
 bool isUnmapped() { return ((b->core.flag & BAM_FUNMAP) != 0); }
 bool isMapped() { return ((b->core.flag & BAM_FUNMAP) == 0); }
 bool isPaired() { return ((b->core.flag & BAM_FPAIRED) != 0); }
 const char* name() { return bam1_qname(b); }
 int pairOrder() {
    //which read in the pair: 0 = unpaired, 1=first read, 2=second read
    int r=0;
    if ((b->core.flag & BAM_FREAD1) != 0) r=1;
    else if ((b->core.flag & BAM_FREAD2) != 0) r=2;
    return r;
    }
 bool revStrand() {
   //this is the raw alignment strand, NOT the transcription (XS) strand
   return ((b->core.flag & BAM_FREVERSE) != 0);
   }
 const char* refName() {
   return (bam_header!=NULL) ?
         ((b->core.tid<0) ? "*" : bam_header->target_name[b->core.tid]) : NULL;
   }
 int32_t refId() { return b->core.tid; }
 int32_t mate_refId() { return b->core.mtid; }
 const char* mate_refName() {
    return (bam_header!=NULL) ?
       ((b->core.mtid<0) ? "*" : bam_header->target_name[b->core.mtid]) : NULL;
    }
 int32_t insertSize() { return b->core.isize; }
 int32_t mate_start() { return b->core.mpos<0? 0 : b->core.mpos+1; }

 int find_tag(const char tag[2], uint8_t* & s, char& tag_type);
 //position s at the beginning of tag data, tag_type is set to the found tag type
 //returns length of tag data, or 0 if tag not found

 char* tag_str(const char tag[2]); //return tag value for tag type 'Z'
 int tag_int(const char tag[2]); //return numeric value of tag (for numeric types)
 char tag_char(const char tag[2]); //return char value of tag (for type 'A')
 char spliceStrand(); // '+', '-' from the XS tag, or '.' if no XS tag

 char* sequence(); //user should free after use
 char* qualities();//user should free after use
 char* cigar(); //returns text version of the CIGAR string; user must free
};


class GBamReader {
   samfile_t* bam_file;
   char* fname;
 public:
   void bopen(const char* filename) {
      if (strcmp(filename, "-")==0) {
        //if stdin was given, we HAVE to assume it's text SAM format, sorry
        bam_file=samopen(filename, "r", 0);
        }
      else {
      //BAM files have the zlib signature bytes at the beginning: 1F 8B 08
      //if that's not present, we assume sam file
      FILE* f=fopen(filename, "rb");
      if (f==NULL) GError("Error opening file %s!\n", filename);
      byte fsig[3];
      size_t rd=fread(fsig, 1, 3, f);
      if (rd<3) GError("Error reading from file %s!\n",filename);
      if ((fsig[0]==0x1F && fsig[1]==0x8B && fsig[2]==0x08) ||
          (fsig[0]=='B' && fsig[1]=='A' && fsig[2]=='M')) {
          bam_file=samopen(filename, "rb", 0); //BAM or uncompressed BAM
          }
        else { //assume text SAM file
          bam_file=samopen(filename, "r", 0);
          }
        }
      if (bam_file==NULL)
         GError("Error: could not open SAM file %s!\n",filename);
      fname=Gstrdup(filename);
      }
   GBamReader(const char* fn) {
      bam_file=NULL;
      fname=NULL;
      bopen(fn);
      }

   bam_header_t* header() {
      return bam_file? bam_file->header : NULL;
      }
   void bclose() {
      if (bam_file) {
        samclose(bam_file);
        bam_file=NULL;
        }
      }
   ~GBamReader() {
      bclose();
      GFREE(fname);
      }
   void rewind() {
     if (fname==NULL) {
       GMessage("Warning: GBamReader::rewind() called without a file name.\n");
       return;
       }
     bclose();
     char* ifname=fname;
     bopen(ifname);
     GFREE(ifname);
     }

   GBamRecord* next() {
      if (bam_file==NULL)
        GError("Warning: GBamReader::next() called with no open file.\n");
      bam1_t* b = bam_init1();
      if (samread(bam_file, b) >= 0) {
        GBamRecord* bamrec=new GBamRecord(b, bam_file->header);
        return bamrec;
        }
      else {
        bam_destroy1(b);
        return NULL;
        }
      }
};


class GBamWriter {
   samfile_t* bam_file;
   bam_header_t* bam_header;
 public:
   void create(const char* fname, bool uncompressed=false) {
      if (bam_header==NULL)
         GError("Error: no bam_header for GBamWriter::create()!\n");
      if (uncompressed) {
         bam_file=samopen(fname, "wbu", bam_header);
         }
        else {
         bam_file=samopen(fname, "wb", bam_header);
         }
      if (bam_file==NULL)
         GError("Error: could not create BAM file %s!\n",fname);
      }
   void create(const char* fname, bam_header_t* bh, bool uncompressed=false) {
     bam_header=bh;
     create(fname,uncompressed);
     }

   GBamWriter(const char* fname, bam_header_t* bh, bool uncompressed=false) {
      create(fname, bh, uncompressed);
      }

   GBamWriter(const char* fname, const char* samfname, bool uncompressed=false) {
      tamFile samf_in=sam_open(samfname);
      if (samf_in==NULL)
         GError("Error: could not open SAM file %s\n", samfname);
      bam_header=sam_header_read(samf_in);
      if (bam_header==NULL)
         GError("Error: could not read SAM header from %s!\n",samfname);
      sam_close(samf_in);
      create(fname, uncompressed);
      }

    ~GBamWriter() {
      samclose(bam_file);
      bam_header_destroy(bam_header);
      }
   bam_header_t* get_header() { return bam_header; }
   int32_t get_tid(const char *seq_name) {
      if (bam_header==NULL)
         GError("Error: missing SAM header (get_tid())\n");
      return bam_get_tid(bam_header, seq_name);
      }

   //just a convenience function for creating a new record, but it's NOT written
   //given pos must be 1-based (so it'll be stored as pos-1 because BAM is 0-based)
   GBamRecord* new_record(const char* qname, const char* gseqname,
            int pos, bool reverse, const char* qseq, const char* cigar=NULL, const char* qual=NULL) {
      int32_t gseq_tid=get_tid(gseqname);
      if (gseq_tid < 0 && strcmp(gseqname, "*")) {
            if (bam_header->n_targets == 0) {
               GError("Error: missing/invalid SAM header\n");
               } else
                   GMessage("Warning: reference '%s' not found in header, will consider it '*'.\n",
                                   gseqname);
            }

      return (new GBamRecord(qname, gseq_tid, pos, reverse, qseq, cigar, qual));
      }

   GBamRecord* new_record(const char* qname, int32_t flags, const char* gseqname,
         int pos, int map_qual, const char* cigar, const char* mgseqname, int mate_pos,
         int insert_size, const char* qseq, const char* quals=NULL,
                          GVec<char*>* aux_strings=NULL) {
      int32_t gseq_tid=get_tid(gseqname);
      if (gseq_tid < 0 && strcmp(gseqname, "*")) {
            if (bam_header->n_targets == 0) {
               GError("Error: missing/invalid SAM header\n");
               } else
                   GMessage("Warning: reference '%s' not found in header, will consider it '*'.\n",
                                   gseqname);
            }
      int32_t mgseq_tid=-1;
      if (mgseqname!=NULL) {
         if (strcmp(mgseqname, "=")==0) {
            mgseq_tid=gseq_tid;
            }
          else {
            mgseq_tid=get_tid(mgseqname);
            if (mgseq_tid < 0 && strcmp(mgseqname, "*")) {
                GMessage("Warning: reference '%s' not found in header, will consider it '*'.\n",
                                   mgseqname);
                }
            }
          }
      return (new GBamRecord(qname, flags, gseq_tid, pos, map_qual, cigar,
              mgseq_tid, mate_pos, insert_size, qseq, quals, aux_strings));
      }

   void write(GBamRecord* brec) {
      if (brec!=NULL)
          samwrite(this->bam_file,brec->get_b());
      }
   void write(bam1_t* b) {
      samwrite(this->bam_file, b);
      }
};

#endif
Revision:	144
Committed:	Thu Dec 22 19:15:46 2011 UTC (12 years, 8 months ago) by gpertea
File size:	11467 byte(s)
Log Message:	added GBitVec.h; removed operator> requirement for sorted vectors; swap() is now template Gswap()
Line	File contents
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