1 static char rcsid[] = "$Id: block.c 218152 2019-01-17 05:37:56Z twu $";
2 #ifdef HAVE_CONFIG_H
3 #include <config.h>
4 #endif
5
6 #include "block.h"
7
8 #include <stdio.h>
9 #include <stdlib.h>
10 #include "mem.h"
11 #ifdef PMAP
12 #include "oligop.h"
13 #else
14 #include "oligo.h"
15 #endif
16 #include "reader.h" /* also has cDNAEnd_T */
17
18 #ifdef DEBUG
19 #define debug(x) x
20 #else
21 #define debug(x)
22 #endif
23
24
25 /* In reading inward from both ends, the block5 and block3 share the
26 same reader. But in reading outward from the middle, the block5
27 and block3 each have their own reader. */
28
29 #define T Block_T
30 struct T {
31 Reader_T reader;
32
33 cDNAEnd_T cdnaend;
34
35 int querylength;
36 int last_querypos;
37 Oligostate_T last_state;
38 Width_T oligosize;
39
40 #ifdef PMAP
41 Oligospace_T aaindex;
42 #else
43 int leftreadshift;
44 Oligospace_T forward;
45 Oligospace_T revcomp;
46 Oligospace_T oligomask;
47 #endif
48
49 /* For saving state */
50 int last_querypos_save;
51 Oligostate_T last_state_save;
52 #ifdef PMAP
53 Oligospace_T aaindex_save;
54 #else
55 Oligospace_T forward_save;
56 Oligospace_T revcomp_save;
57 #endif
58 };
59
60
61 Reader_T
Block_reader(T this)62 Block_reader (T this) {
63 return this->reader;
64 }
65
66 int
Block_querypos(T this)67 Block_querypos (T this) {
68 return this->last_querypos;
69 }
70
71 #ifdef PMAP
72 Oligospace_T
Block_aaindex(T this)73 Block_aaindex (T this) {
74 return this->aaindex;
75 }
76 #else
77 Oligospace_T
Block_forward(T this)78 Block_forward (T this) {
79 if (this->cdnaend == FIVE) {
80 return this->forward & this->oligomask;
81 } else {
82 return (this->forward >> this->leftreadshift) & this->oligomask;
83 }
84 }
85
86 Oligospace_T
Block_revcomp(T this)87 Block_revcomp (T this) {
88 if (this->cdnaend == FIVE) {
89 return (this->revcomp >> this->leftreadshift) & this->oligomask;
90 } else {
91 return this->revcomp & this->oligomask;
92 }
93 }
94 #endif
95
96 bool
Block_donep(T this)97 Block_donep (T this) {
98 if (this->last_state == DONE) {
99 return true;
100 } else {
101 return false;
102 }
103 }
104
105
106 extern void
Block_save(T this)107 Block_save (T this) {
108 /* save->reader = this->reader; -- not necessary */
109 /* save->cdnaend = this->cdnaend; -- not necessary */
110
111 this->last_querypos_save = this->last_querypos;
112 this->last_state_save = this->last_state;
113 #ifdef PMAP
114 this->aaindex_save = this->aaindex;
115 debug(printf("Saving block at last_querypos %d, aaindex %u\n",this->last_querypos,this->aaindex));
116 #else
117 this->forward_save = this->forward;
118 this->revcomp_save = this->revcomp;
119 debug(printf("Saving block at last_querypos %d, forward %u, revcomp %u\n",
120 this->last_querypos,this->forward,this->revcomp));
121 #endif
122
123 return;
124 }
125
126 extern void
Block_restore(T this)127 Block_restore (T this) {
128 /* this->reader = save->reader; -- not necessary */
129 /* this->cdnaend = save->cdnaend; -- not necessary */
130
131 if (this->cdnaend == FIVE) {
132 Reader_reset_start(this->reader,this->last_querypos_save+this->oligosize);
133 } else {
134 Reader_reset_end(this->reader,this->last_querypos_save-1);
135 }
136
137 this->last_querypos = this->last_querypos_save;
138 this->last_state = this->last_state_save;
139 #ifdef PMAP
140 this->aaindex = this->aaindex_save;
141 #else
142 this->forward = this->forward_save;
143 this->revcomp = this->revcomp_save;
144 #endif
145 return;
146 }
147
148 extern void
Block_reset_ends(T this)149 Block_reset_ends (T this) {
150 Reader_reset_ends(this->reader);
151
152 if (this->cdnaend == FIVE) {
153 this->last_querypos = Reader_querystart(this->reader);
154 } else {
155 this->last_querypos = Reader_queryend(this->reader);
156 }
157 this->last_state = INIT;
158
159 #ifdef PMAP
160 this->aaindex = 0U;
161 #else
162 this->forward = 0U;
163 this->revcomp = 0U;
164 #endif
165
166 return;
167 }
168
169
170 T
Block_new(cDNAEnd_T cdnaend,Width_T oligosize,int leftreadshift,Reader_T reader,int querylength)171 Block_new (cDNAEnd_T cdnaend, Width_T oligosize,
172 #ifndef PMAP
173 int leftreadshift,
174 #endif
175 Reader_T reader, int querylength) {
176 T new = (T) MALLOC(sizeof(*new));
177
178 new->reader = reader;
179 new->cdnaend = cdnaend;
180 new->oligosize = oligosize;
181 #ifndef PMAP
182 #ifdef HAVE_64_BIT
183 new->oligomask = ~(~0ULL << 2*oligosize);
184 #else
185 new->oligomask = ~(~0U << 2*oligosize);
186 #endif
187 new->leftreadshift = leftreadshift;
188 #endif
189
190 new->querylength = querylength;
191 if (cdnaend == FIVE) {
192 #ifdef PMAP
193 new->last_querypos = Reader_startpos(reader);
194 #else
195 new->last_querypos = Reader_startpos(reader) - oligosize;
196 #endif
197 } else if (cdnaend == THREE) {
198 new->last_querypos = Reader_endpos(reader) + 1;
199 }
200 new->last_state = INIT;
201
202 #ifdef PMAP
203 new->aaindex = 0U;
204 #else
205 new->forward = 0U;
206 new->revcomp = 0U;
207 #endif
208
209 new->last_querypos_save = new->last_querypos;
210 new->last_state_save = new->last_state;
211 #ifdef PMAP
212 new->aaindex_save = new->aaindex;
213 #else
214 new->forward_save = new->forward;
215 new->revcomp_save = new->revcomp;
216 #endif
217
218 return new;
219 }
220
221 void
Block_free(T * old)222 Block_free (T *old) {
223 if (*old) {
224 FREE(*old);
225 }
226 return;
227 }
228
229
230 bool
Block_next_5(T this,int genestrand)231 Block_next_5 (T this, int genestrand) {
232 #ifdef DEBUG
233 char *nt_fwd, *nt_rev;
234 #endif
235
236 if (this->last_state == DONE) {
237 return false;
238 } else {
239
240 #ifdef PMAP
241 this->last_state = Oligo_next_5(this->last_state,&this->last_querypos,
242 &this->aaindex,this->reader,genestrand);
243 debug(printf("Block has aaindex %u at querypos %d\n",
244 this->aaindex,this->last_querypos));
245 #else
246 this->last_state = Oligo_next_5(this->last_state,&this->last_querypos,
247 &this->forward,&this->revcomp,this->reader,genestrand);
248 debug(
249 if (this->cdnaend == THREE) {
250 nt_fwd = Oligo_one_nt(this->forward >> this->leftreadshift,this->oligosize);
251 nt_rev = Oligo_one_nt(this->revcomp,this->oligosize);
252 }
253 if (this->cdnaend == FIVE) {
254 nt_fwd = Oligo_one_nt(this->forward,this->oligosize);
255 nt_rev = Oligo_one_nt(this->revcomp >> this->leftreadshift,this->oligosize);
256 }
257 printf("Block has oligo forward %s, revcomp %s at querypos %d\n",
258 nt_fwd,nt_rev,this->last_querypos);
259 FREE(nt_rev);
260 FREE(nt_fwd)
261 );
262 #endif
263
264 if (this->last_state == DONE) {
265 return false;
266 } else {
267 return true;
268 }
269 }
270 }
271
272
273 bool
Block_next_3(T this,int genestrand)274 Block_next_3 (T this, int genestrand) {
275 #ifdef DEBUG
276 char *nt_fwd, *nt_rev;
277 #endif
278
279 if (this->last_state == DONE) {
280 return false;
281 } else {
282
283 #ifdef PMAP
284 this->last_state = Oligo_next_3(this->last_state,&this->last_querypos,
285 &this->aaindex,this->reader,genestrand);
286 debug(printf("Block has aaindex %u at querypos %d\n",
287 this->aaindex,this->last_querypos));
288 #else
289 this->last_state = Oligo_next_3(this->last_state,&this->last_querypos,
290 &this->forward,&this->revcomp,this->reader,genestrand);
291 debug(
292 if (this->cdnaend == THREE) {
293 nt_fwd = Oligo_one_nt(this->forward >> this->leftreadshift,this->oligosize);
294 nt_rev = Oligo_one_nt(this->revcomp,this->oligosize);
295 }
296 if (this->cdnaend == FIVE) {
297 nt_fwd = Oligo_one_nt(this->forward,this->oligosize);
298 nt_rev = Oligo_one_nt(this->revcomp >> this->leftreadshift,this->oligosize);
299 }
300 printf("Block has oligo forward %s, revcomp %s at querypos %d\n",
301 nt_fwd,nt_rev,this->last_querypos);
302 FREE(nt_rev);
303 FREE(nt_fwd)
304 );
305 #endif
306
307 if (this->last_state == DONE) {
308 return false;
309 } else {
310 return true;
311 }
312 }
313 }
314
315
316 /* Returns whether skipping was successful (as opposed to Block_next, which returns whether scanning can continue) */
317 bool
Block_skip_5(T this,int genestrand,int nskip)318 Block_skip_5 (T this, int genestrand, int nskip) {
319 int init_querypos;
320
321 if (this->last_state == DONE) {
322 return false;
323 } else {
324
325 init_querypos = this->last_querypos;
326 #ifdef PMAP
327 this->last_state = Oligo_skip_5(this->last_state,&this->last_querypos,
328 &this->aaindex,this->reader,genestrand,nskip);
329 debug(printf("Block has aaindex %u at querypos %d\n",
330 this->aaindex,this->last_querypos));
331 #else
332 this->last_state = Oligo_skip_5(this->last_state,&this->last_querypos,
333 &this->forward,&this->revcomp,
334 this->reader,genestrand,nskip);
335 debug(printf("Block has oligo %08X at querypos %d\n",
336 this->forward,this->last_querypos));
337 #endif
338
339 if (this->last_state != VALID) {
340 return false;
341 } else if (this->last_querypos - init_querypos != nskip) {
342 return false;
343 } else {
344 return true;
345 }
346 }
347 }
348
349
350 bool
Block_skip_3(T this,int genestrand,int nskip)351 Block_skip_3 (T this, int genestrand, int nskip) {
352 int init_querypos;
353
354 if (this->last_state == DONE) {
355 return false;
356 } else {
357
358 init_querypos = this->last_querypos;
359 #ifdef PMAP
360 this->last_state = Oligo_skip_3(this->last_state,&this->last_querypos,
361 &this->aaindex,this->reader,genestrand,nskip);
362 debug(printf("Block has aaindex %u at querypos %d\n",
363 this->aaindex,this->last_querypos));
364 #else
365 this->last_state = Oligo_skip_3(this->last_state,&this->last_querypos,
366 &this->forward,&this->revcomp,
367 this->reader,genestrand,nskip);
368 debug(printf("Block has oligo %08X at querypos %d\n",
369 this->forward,this->last_querypos));
370 #endif
371
372 if (this->last_state != VALID) {
373 return false;
374 } else if (this->last_querypos - init_querypos != nskip) {
375 return false;
376 } else {
377 return true;
378 }
379 }
380 }
381
382
383 #if 0
384 /* Returns whether skipping was successful (as opposed to Block_next, which returns whether scanning can continue) */
385 bool
386 Block_skipto (T this, int querypos) {
387 int nskip;
388
389 if (this->cdnaend == FIVE) {
390 nskip = querypos - this->last_querypos;
391 } else {
392 nskip = this->last_querypos - querypos;
393 }
394 if (nskip < 0) {
395 return false;
396 } else {
397 return Block_skip(this,nskip);
398 }
399 }
400 #endif
401
402
403 /* Returns querypos */
404 #ifdef PMAP
405
406 int
Block_process_oligo(Univcoord_T ** fwdpositions,int * nfwdhits,Univcoord_T ** revpositions,int * nrevhits,T this,Indexdb_T indexdb_fwd,Indexdb_T indexdb_rev,int indexdb_sizelimit)407 Block_process_oligo (Univcoord_T **fwdpositions, int *nfwdhits,
408 Univcoord_T **revpositions, int *nrevhits,
409 T this, Indexdb_T indexdb_fwd, Indexdb_T indexdb_rev,
410 int indexdb_sizelimit) {
411
412 /* Note that querylength was already multiplied by 3 */
413 *fwdpositions = Indexdb_read_with_diagterm(&(*nfwdhits),indexdb_fwd,this->aaindex,
414 /*diagterm*/this->querylength-3*this->last_querypos);
415 *revpositions = Indexdb_read_with_diagterm(&(*nrevhits),indexdb_rev,this->aaindex,
416 /*diagterm*/3*this->last_querypos);
417 return this->last_querypos;
418 }
419
420 #else
421
422 /* In standard mode, indexdb_rev is indexdb_fwd. They differ for cmet and atoi modes. */
423 int
Block_process_oligo_5(Univcoord_T ** fwdpositions,int * nfwdhits,Univcoord_T ** revpositions,int * nrevhits,T this,Indexdb_T indexdb_fwd,Indexdb_T indexdb_rev)424 Block_process_oligo_5 (Univcoord_T **fwdpositions, int *nfwdhits,
425 Univcoord_T **revpositions, int *nrevhits,
426 T this, Indexdb_T indexdb_fwd, Indexdb_T indexdb_rev) {
427 #if 0
428 printf("block_process: Querypos %d, oligos are %06X and %06X\n",this->last_querypos,this->forward,this->revcomp >> this->leftreadshift);
429 #endif
430 *fwdpositions = Indexdb_read_with_diagterm(&(*nfwdhits),indexdb_fwd,this->forward & this->oligomask,
431 /*diagterm*/this->querylength-this->last_querypos);
432 *revpositions = Indexdb_read_with_diagterm(&(*nrevhits),indexdb_rev,(this->revcomp >> this->leftreadshift) & this->oligomask,
433 /*diagterm*/this->last_querypos);
434
435 return this->last_querypos;
436 }
437
438 /* In standard mode, indexdb_rev is indexdb_fwd. They differ for cmet and atoi modes. */
439 int
Block_process_oligo_3(Univcoord_T ** fwdpositions,int * nfwdhits,Univcoord_T ** revpositions,int * nrevhits,T this,Indexdb_T indexdb_fwd,Indexdb_T indexdb_rev)440 Block_process_oligo_3 (Univcoord_T **fwdpositions, int *nfwdhits,
441 Univcoord_T **revpositions, int *nrevhits,
442 T this, Indexdb_T indexdb_fwd, Indexdb_T indexdb_rev) {
443 #if 0
444 printf("block_process Querypos %d, oligos are %06X and %06X\n",this->last_querypos,this->forward >> this->leftreadshift,this->revcomp);
445 #endif
446 *fwdpositions = Indexdb_read_with_diagterm(&(*nfwdhits),indexdb_fwd,(this->forward >> this->leftreadshift) & this->oligomask,
447 /*diagterm*/this->querylength-this->last_querypos);
448 *revpositions = Indexdb_read_with_diagterm(&(*nrevhits),indexdb_rev,this->revcomp & this->oligomask,
449 /*diagterm*/this->last_querypos);
450
451 return this->last_querypos;
452 }
453
454
455 #endif
456