annotate rDiff/mex/get_reads_direct.cpp @ 0:0f80a5141704

version 0.3 uploaded
author vipints
date Thu, 14 Feb 2013 23:38:36 -0500
parents
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
0
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
1 /* written by Jonas Behr, Regina Bohnert and Gunnar Raetsch, FML Tuebingen, Germany, 2010 */
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
2
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
3 #include <stdio.h>
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
4 #include <assert.h>
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
5 #include "sam.h"
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
6 #include "get_reads_direct.h"
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
7
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
8 #include <vector>
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
9 using std::vector;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
10 #include <string>
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
11 using std::string;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
12
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
13 #ifdef __cplusplus
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
14 extern "C" {
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
15 #endif
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
16
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
17 typedef struct {
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
18 int beg, end;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
19 samfile_t *in;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
20 } tmpstruct_t;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
21
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
22 typedef struct {
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
23 uint64_t u, v;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
24 } pair64_t;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
25
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
26 static inline int is_overlap(uint32_t beg, uint32_t end, const bam1_t *b)
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
27 {
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
28 uint32_t rbeg = b->core.pos;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
29 uint32_t rend = b->core.n_cigar? bam_calend(&b->core, bam1_cigar(b)) : b->core.pos + 1;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
30 return (rend > beg && rbeg < end);
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
31 }
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
32
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
33 pair64_t * get_chunk_coordinates(const bam_index_t *idx, int tid, int beg, int end, int* cnt_off);
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
34
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
35 int bam_fetch_reads(bamFile fp, const bam_index_t *idx, int tid, int beg, int end, void *data, bam_header_t* header, vector<CRead*>* reads, char strand);
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
36
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
37 // callback for bam_plbuf_init()
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
38 static int pileup_func(uint32_t tid, uint32_t pos, int n, const bam_pileup1_t *pl, void *data)
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
39 {
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
40 //tmpstruct_t *tmp = (tmpstruct_t*)data;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
41 //if ((int)pos >= tmp->beg && (int)pos < tmp->end)
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
42 // printf("%s\t%d\t%d\n", tmp->in->header->target_name[tid], pos + 1, n);
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
43 return 0;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
44 }
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
45 #ifdef __cplusplus
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
46 }
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
47 #endif
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
48 int parse_sam_line(char* line, CRead* read);
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
49 int set_strand(char c);
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
50 //void parse_cigar(bam1_t* b, CRead* read);
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
51 void parse_cigar(bam1_t* b, CRead* read, bam_header_t* header);
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
52
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
53
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
54 int get_reads_from_bam(char* filename, char* region, vector<CRead*>* reads, char strand, int lsubsample)
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
55 {
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
56 subsample = lsubsample;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
57 set_strand(strand);
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
58
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
59 srand (time(NULL));
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
60 //srand (1234);
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
61 tmpstruct_t tmp;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
62 tmp.in = samopen(filename, "rb", 0);
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
63 if (tmp.in == 0) {
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
64 fprintf(stderr, "Fail to open BAM file %s\n", filename);
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
65 return 1;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
66 }
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
67 int ref;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
68 bam_index_t *idx;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
69 bam_plbuf_t *buf;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
70 idx = bam_index_load(filename); // load BAM index
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
71 if (idx == 0) {
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
72 fprintf(stderr, "BAM indexing file is not available.\n");
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
73 return 1;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
74 }
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
75 bam_parse_region(tmp.in->header, region, &ref,
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
76 &tmp.beg, &tmp.end); // parse the region
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
77 if (ref < 0) {
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
78 fprintf(stderr, "Invalid region %s\n", region);
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
79 return 1;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
80 }
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
81
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
82 buf = bam_plbuf_init(pileup_func, &tmp); // initialize pileup
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
83
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
84 bam_fetch_reads(tmp.in->x.bam, idx, ref, tmp.beg, tmp.end, buf, tmp.in->header, reads, strand);
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
85 //fprintf(stdout, "intron_list: %d \n", intron_list->size());
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
86
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
87 bam_plbuf_push(0, buf); // finalize pileup
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
88 bam_index_destroy(idx);
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
89 bam_plbuf_destroy(buf);
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
90 samclose(tmp.in);
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
91 return 0;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
92 }
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
93
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
94
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
95 int bam_fetch_reads(bamFile fp, const bam_index_t *idx, int tid, int beg, int end, void *data, bam_header_t* header, vector<CRead*>* reads, char strand)
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
96 {
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
97 int n_off;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
98 pair64_t *off = get_chunk_coordinates(idx, tid, beg, end, &n_off);
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
99 if (off == 0) return 0;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
100 {
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
101 // retrive alignments
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
102 uint64_t curr_off;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
103 int i, ret, n_seeks;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
104 n_seeks = 0; i = -1; curr_off = 0;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
105 bam1_t *b = (bam1_t*)calloc(1, sizeof(bam1_t));
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
106 for (;;) {
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
107 if (curr_off == 0 || curr_off >= off[i].v) { // then jump to the next chunk
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
108 if (i == n_off - 1) break; // no more chunks
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
109 if (i >= 0) assert(curr_off == off[i].v); // otherwise bug
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
110 if (i < 0 || off[i].v != off[i+1].u) { // not adjacent chunks; then seek
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
111 bam_seek(fp, off[i+1].u, SEEK_SET);
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
112 curr_off = bam_tell(fp);
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
113 ++n_seeks;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
114 }
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
115 ++i;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
116 }
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
117 if ((ret = bam_read1(fp, b)) > 0) {
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
118 curr_off = bam_tell(fp);
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
119 if (b->core.tid != tid || b->core.pos >= end) break; // no need to proceed
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
120 else if (is_overlap(beg, end, b))
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
121 {
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
122 int rr = rand();
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
123 if ((rr%1000 < subsample))
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
124 {
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
125 CRead* read = new CRead();
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
126 parse_cigar(b, read, header);
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
127
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
128 if (strand == '0' || strand==read->strand[0])
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
129 {
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
130 reads->push_back(read);
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
131 }
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
132 else if (read->strand[0]=='0')
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
133 {
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
134 reads->push_back(read);
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
135 }
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
136 //else if (read->strand[0]=='0'&&((b->core.flag & g_flag_off) >0))
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
137 //{
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
138 // //fprintf(stdout, "(-)-strand; read->strand[0]==0, num_exons: %i \n", read->block_starts.size());
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
139 // // this flag means that the read has been reversed for alignment
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
140 // // flag bit set and (-)-strand requested
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
141 // reads->push_back(read);
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
142 //}
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
143 //else if (read->strand[0]=='0'&&(g_flag_on>0&&(b->core.flag & g_flag_on)==0))
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
144 //{
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
145 // //fprintf(stdout, "(+)-strand; read->strand[0]==0, num_exons: %i \n", read->block_starts.size());
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
146 // // (+)-strand requested and flag bit not set
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
147 // reads->push_back(read);
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
148 //}
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
149 }
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
150 }
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
151 } else break; // end of file
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
152 }
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
153 // fprintf(stderr, "[bam_fetch] # seek calls: %d\n", n_seeks);
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
154 bam_destroy1(b);
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
155 }
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
156 free(off);
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
157 return 0;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
158 }
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
159
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
160 void parse_cigar(bam1_t* b, CRead* read, bam_header_t* header)
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
161 {
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
162 read->start_pos = b->core.pos+1;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
163 read->set_strand('0');
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
164
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
165 for (int k = 0; k < b->core.n_cigar; ++k)
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
166 {
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
167 int op = bam1_cigar(b)[k] & BAM_CIGAR_MASK; // operation
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
168 int l = bam1_cigar(b)[k] >> BAM_CIGAR_SHIFT; // length
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
169 //fprintf(stdout, "op:%d l:%d\n", op, l);
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
170 if (op == BAM_CMATCH)
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
171 {
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
172 if (k==0)
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
173 {
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
174 read->block_lengths.push_back(l);
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
175 read->block_starts.push_back(0);
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
176 }
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
177 else
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
178 {
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
179 int op_prev = bam1_cigar(b)[k-1] & BAM_CIGAR_MASK;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
180 int l_prev = bam1_cigar(b)[k-1] >> BAM_CIGAR_SHIFT;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
181 if (op_prev==BAM_CREF_SKIP)// intron before
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
182 {
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
183 if (read->block_lengths.size()>=1)
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
184 {
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
185 int last_block_start = (*(read->block_starts.end()-1));
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
186 int intron_start = last_block_start+(*(read->block_lengths.end()-1));
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
187 read->block_lengths.push_back(l);
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
188 read->block_starts.push_back(intron_start+l_prev);
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
189 }
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
190 else
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
191 {
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
192 // start of first block was not a match
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
193 read->block_lengths.push_back(l);
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
194 read->block_starts.push_back(0);
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
195 }
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
196 }
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
197 else
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
198 {
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
199 if (read->block_lengths.size()>=1 && op == BAM_CDEL)// if it is an insertion then the matching block is not inreased
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
200 (*(read->block_lengths.end()-1))+=l;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
201 else
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
202 {
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
203 //char *samline = bam_format1(header, b);
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
204 //printf("header: %s \n", samline);
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
205 }
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
206 }
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
207 }
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
208 }
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
209 else if (op == BAM_CDEL)
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
210 {
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
211 if (k>0 && read->block_lengths.size()>=1)
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
212 (*(read->block_lengths.end()-1))+=l;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
213 }
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
214 else if (op == BAM_CREF_SKIP)//intron
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
215 {}
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
216 else if (op == BAM_CINS || op == BAM_CSOFT_CLIP)
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
217 {}
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
218 }
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
219 // parse auxiliary data
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
220 uint8_t* s = bam1_aux(b);
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
221 uint8_t* end = b->data + b->data_len;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
222 while (s < end)
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
223 {
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
224 uint8_t type, key[2];
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
225 key[0] = s[0]; key[1] = s[1];
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
226 s += 2; type = *s; ++s;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
227 //fprintf(stdout, "\n%c%c:%c\n", key[0], key[1], type);
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
228 if (type == 'A')
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
229 {
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
230 if ( key[0] =='X' && key[1] == 'S')
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
231 {
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
232 read->set_strand((char) *s);
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
233 }
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
234 ++s;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
235 }
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
236 else if (type == 'C')
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
237 {
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
238 if ( key[0] =='H' && key[1] == '0')
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
239 {
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
240 uint8_t matches = *s;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
241 read->matches = (int) matches;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
242 }
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
243 if ( key[0] =='N' && key[1] == 'M')
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
244 {
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
245 uint8_t mismatches = *s;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
246 read->mismatches = (int) mismatches;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
247 }
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
248 if ( key[0] =='H' && key[1] == 'I')
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
249 {
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
250 uint8_t mai = *s;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
251 read->multiple_alignment_index = (int) mai;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
252 }
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
253
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
254 ++s;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
255 }
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
256 else if (type == 'c') { ++s; }
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
257 else if (type == 'S') { s += 2; }
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
258 else if (type == 's') { s += 2; }
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
259 else if (type == 'I') { s += 4; }
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
260 else if (type == 'i') { s += 4; }
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
261 else if (type == 'f') { s += 4; }
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
262 else if (type == 'd') { s += 8; }
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
263 else if (type == 'Z') { ++s; }
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
264 else if (type == 'H') { ++s; }
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
265 }
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
266 }
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
267
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
268 int set_strand(char c)
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
269 {
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
270 if (c=='+')
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
271 {
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
272 char* fl = (char*) "0x0010";
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
273 g_flag_on = strtol(fl, 0, 0);
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
274 g_flag_off = 0;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
275 }
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
276 else if (c=='-')
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
277 {
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
278 char* fl = (char*) "0x0010";
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
279 g_flag_off = strtol(fl, 0, 0);
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
280 g_flag_on = 0;
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
281 }
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
282 }
0f80a5141704 version 0.3 uploaded
vipints
parents:
diff changeset
283