annotate pyPRADA_1.2/tools/bwa-0.5.7-mh/bwtsw2_aux.c @ 0:acc2ca1a3ba4

Uploaded
author siyuan
date Thu, 20 Feb 2014 00:44:58 -0500
parents
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
0
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
1 #include <stdlib.h>
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
2 #include <stdio.h>
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
3 #include <math.h>
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
4 #ifdef HAVE_CONFIG_H
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
5 #include "config.h"
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
6 #endif
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
7 #ifdef HAVE_PTHREAD
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
8 #include <pthread.h>
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
9 #endif
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
10 #include "bntseq.h"
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
11 #include "bwt_lite.h"
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
12 #include "utils.h"
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
13 #include "bwtsw2.h"
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
14 #include "stdaln.h"
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
15 #include "kstring.h"
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
16
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
17 #include "kseq.h"
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
18 KSEQ_INIT(gzFile, gzread)
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
19
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
20 #include "ksort.h"
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
21 #define __left_lt(a, b) ((a).end > (b).end)
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
22 KSORT_INIT(hit, bsw2hit_t, __left_lt)
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
23
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
24 extern unsigned char nst_nt4_table[256];
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
25
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
26 unsigned char nt_comp_table[256] = {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
27 'N','N','N','N', 'N','N','N','N', 'N','N','N','N', 'N','N','N','N',
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
28 'N','N','N','N', 'N','N','N','N', 'N','N','N','N', 'N','N','N','N',
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
29 'N','N','N','N', 'N','N','N','N', 'N','N','N','N', 'N','N','N','N',
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
30 'N','N','N','N', 'N','N','N','N', 'N','N','N','N', 'N','N','N','N',
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
31 'N','T','V','G', 'H','N','N','C', 'D','N','N','M', 'N','K','N','N',
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
32 'N','N','Y','S', 'A','N','B','W', 'X','R','N','N', 'N','N','N','N',
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
33 'n','t','v','g', 'h','n','n','c', 'd','n','n','m', 'n','k','n','n',
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
34 'n','n','y','s', 'a','n','b','w', 'x','r','n','N', 'N','N','N','N',
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
35 'N','N','N','N', 'N','N','N','N', 'N','N','N','N', 'N','N','N','N',
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
36 'N','N','N','N', 'N','N','N','N', 'N','N','N','N', 'N','N','N','N',
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
37 'N','N','N','N', 'N','N','N','N', 'N','N','N','N', 'N','N','N','N',
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
38 'N','N','N','N', 'N','N','N','N', 'N','N','N','N', 'N','N','N','N',
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
39 'N','N','N','N', 'N','N','N','N', 'N','N','N','N', 'N','N','N','N',
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
40 'N','N','N','N', 'N','N','N','N', 'N','N','N','N', 'N','N','N','N',
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
41 'N','N','N','N', 'N','N','N','N', 'N','N','N','N', 'N','N','N','N',
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
42 'N','N','N','N', 'N','N','N','N', 'N','N','N','N', 'N','N','N','N'
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
43 };
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
44
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
45 extern int bsw2_resolve_duphits(const bwt_t *bwt, bwtsw2_t *b, int IS);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
46 extern int bsw2_resolve_query_overlaps(bwtsw2_t *b, float mask_level);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
47
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
48 bsw2opt_t *bsw2_init_opt()
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
49 {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
50 bsw2opt_t *o = (bsw2opt_t*)calloc(1, sizeof(bsw2opt_t));
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
51 o->a = 1; o->b = 3; o->q = 5; o->r = 2; o->t = 30;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
52 o->bw = 50;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
53 o->z = 1; o->is = 3; o->t_seeds = 5; o->hard_clip = 0;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
54 o->mask_level = 0.50f; o->yita = 5.5f; o->coef = 5.5f;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
55 o->qr = o->q + o->r; o->n_threads = 1; o->chunk_size = 10000000;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
56 return o;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
57 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
58
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
59 void bsw2_destroy(bwtsw2_t *b)
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
60 {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
61 int i;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
62 if (b == 0) return;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
63 if (b->cigar)
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
64 for (i = 0; i < b->n; ++i) free(b->cigar[i]);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
65 free(b->cigar); free(b->n_cigar); free(b->hits);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
66 free(b);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
67 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
68
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
69 #define __gen_ap(par, opt) do { \
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
70 int i; \
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
71 for (i = 0; i < 25; ++i) (par).matrix[i] = -(opt)->b; \
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
72 for (i = 0; i < 4; ++i) (par).matrix[i*5+i] = (opt)->a; \
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
73 (par).gap_open = (opt)->q; (par).gap_ext = (opt)->r; \
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
74 (par).gap_end = (opt)->r; \
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
75 (par).row = 5; (par).band_width = opt->bw; \
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
76 } while (0)
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
77
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
78 #define __rpac(pac, l, i) (pac[(l-i-1)>>2] >> (~(l-i-1)&3)*2 & 0x3)
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
79
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
80 void bsw2_extend_left(const bsw2opt_t *opt, bwtsw2_t *b, uint8_t *_query, int lq, uint8_t *pac, uint32_t l_pac, int is_rev, uint8_t *_mem)
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
81 {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
82 int i, matrix[25];
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
83 bwtint_t k;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
84 uint8_t *target = 0, *query;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
85 AlnParam par;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
86
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
87 par.matrix = matrix;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
88 __gen_ap(par, opt);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
89 query = calloc(lq, 1);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
90 // sort according to the descending order of query end
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
91 ks_introsort(hit, b->n, b->hits);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
92 target = calloc(((lq + 1) / 2 * opt->a + opt->r) / opt->r + lq, 1);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
93 // reverse _query
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
94 for (i = 0; i < lq; ++i) query[lq - i - 1] = _query[i];
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
95 // core loop
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
96 for (i = 0; i < b->n; ++i) {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
97 bsw2hit_t *p = b->hits + i;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
98 int lt = ((p->beg + 1) / 2 * opt->a + opt->r) / opt->r + lq;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
99 int score, j;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
100 path_t path;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
101 p->n_seeds = 1;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
102 if (p->l || p->k == 0) continue;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
103 for (j = score = 0; j < i; ++j) {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
104 bsw2hit_t *q = b->hits + j;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
105 if (q->beg <= p->beg && q->k <= p->k && q->k + q->len >= p->k + p->len) {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
106 if (q->n_seeds < (1<<14) - 2) ++q->n_seeds;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
107 ++score;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
108 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
109 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
110 if (score) continue;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
111 if (lt > p->k) lt = p->k;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
112 if (is_rev) {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
113 for (k = p->k - 1, j = 0; k > 0 && j < lt; --k) // FIXME: k=0 not considered!
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
114 target[j++] = __rpac(pac, l_pac, k);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
115 } else {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
116 for (k = p->k - 1, j = 0; k > 0 && j < lt; --k) // FIXME: k=0 not considered!
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
117 target[j++] = pac[k>>2] >> (~k&3)*2 & 0x3;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
118 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
119 lt = j;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
120 score = aln_extend_core(target, lt, query + lq - p->beg, p->beg, &par, &path, 0, p->G, _mem);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
121 if (score > p->G) { // extensible
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
122 p->G = score;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
123 p->len += path.i;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
124 p->beg -= path.j;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
125 p->k -= path.i;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
126 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
127 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
128 free(query); free(target);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
129 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
130
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
131 void bsw2_extend_rght(const bsw2opt_t *opt, bwtsw2_t *b, uint8_t *query, int lq, uint8_t *pac, uint32_t l_pac, int is_rev, uint8_t *_mem)
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
132 {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
133 int i, matrix[25];
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
134 uint32_t k;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
135 uint8_t *target;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
136 AlnParam par;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
137
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
138 par.matrix = matrix;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
139 __gen_ap(par, opt);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
140 target = calloc(((lq + 1) / 2 * opt->a + opt->r) / opt->r + lq, 1);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
141 for (i = 0; i < b->n; ++i) {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
142 bsw2hit_t *p = b->hits + i;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
143 int lt = ((lq - p->beg + 1) / 2 * opt->a + opt->r) / opt->r + lq;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
144 int j, score;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
145 path_t path;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
146 if (p->l) continue;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
147 if (is_rev) {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
148 for (k = p->k, j = 0; k < p->k + lt && k < l_pac; ++k)
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
149 target[j++] = __rpac(pac, l_pac, k);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
150 } else {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
151 for (k = p->k, j = 0; k < p->k + lt && k < l_pac; ++k)
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
152 target[j++] = pac[k>>2] >> (~k&3)*2 & 0x3;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
153 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
154 lt = j;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
155 score = aln_extend_core(target, lt, query + p->beg, lq - p->beg, &par, &path, 0, 1, _mem);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
156 // if (score < p->G) fprintf(stderr, "[bsw2_extend_hits] %d < %d\n", score, p->G);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
157 if (score >= p->G) {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
158 p->G = score;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
159 p->len = path.i;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
160 p->end = path.j + p->beg;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
161 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
162 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
163 free(target);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
164 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
165
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
166 /* generate CIGAR array(s) in b->cigar[] */
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
167 static void gen_cigar(const bsw2opt_t *opt, int lq, uint8_t *seq[2], uint8_t *pac, bwtsw2_t *b)
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
168 {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
169 uint8_t *target;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
170 int i, matrix[25];
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
171 AlnParam par;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
172 path_t *path;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
173
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
174 par.matrix = matrix;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
175 __gen_ap(par, opt);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
176 i = ((lq + 1) / 2 * opt->a + opt->r) / opt->r + lq; // maximum possible target length
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
177 target = calloc(i, 1);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
178 path = calloc(i + lq, sizeof(path_t));
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
179 // memory clean up for b
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
180 if (b->n < b->max) {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
181 b->max = b->n;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
182 b->hits = realloc(b->hits, b->n * sizeof(bsw2hit_t));
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
183 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
184 if (b->cigar) free(b->cigar);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
185 if (b->n_cigar) free(b->n_cigar);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
186 b->cigar = (uint32_t**)calloc(b->max, sizeof(void*));
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
187 b->n_cigar = (int*)calloc(b->max, sizeof(int));
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
188 // generate CIGAR
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
189 for (i = 0; i < b->n; ++i) {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
190 bsw2hit_t *p = b->hits + i;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
191 uint8_t *query;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
192 uint32_t k;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
193 int score, path_len, beg, end;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
194 if (p->l) continue;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
195 beg = (p->flag & 0x10)? lq - p->end : p->beg;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
196 end = (p->flag & 0x10)? lq - p->beg : p->end;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
197 query = seq[(p->flag & 0x10)? 1 : 0] + beg;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
198 for (k = p->k; k < p->k + p->len; ++k) // in principle, no out-of-boundary here
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
199 target[k - p->k] = pac[k>>2] >> (~k&3)*2 & 0x3;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
200 score = aln_global_core(target, p->len, query, end - beg, &par, path, &path_len);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
201 b->cigar[i] = aln_path2cigar32(path, path_len, &b->n_cigar[i]);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
202 if (beg != 0 || end < lq) { // write soft clipping
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
203 b->cigar[i] = realloc(b->cigar[i], 4 * (b->n_cigar[i] + 2));
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
204 if (beg != 0) {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
205 memmove(b->cigar[i] + 1, b->cigar[i], b->n_cigar[i] * 4);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
206 b->cigar[i][0] = beg<<4 | 4;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
207 ++b->n_cigar[i];
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
208 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
209 if (end < lq) {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
210 b->cigar[i][b->n_cigar[i]] = (lq - end)<<4 | 4;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
211 ++b->n_cigar[i];
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
212 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
213 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
214 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
215 free(target); free(path);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
216 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
217
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
218 /* this is for the debugging purpose only */
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
219 void bsw2_debug_hits(const bwtsw2_t *b)
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
220 {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
221 int i;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
222 printf("# raw hits: %d\n", b->n);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
223 for (i = 0; i < b->n; ++i) {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
224 bsw2hit_t *p = b->hits + i;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
225 if (p->l == 0)
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
226 printf("%d, %d, %d, %u, %u\n", p->G, p->beg, p->end, p->k, p->l);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
227 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
228 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
229
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
230 static void merge_hits(bwtsw2_t *b[2], int l, int is_reverse)
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
231 {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
232 int i;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
233 if (b[0]->n + b[1]->n > b[0]->max) {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
234 b[0]->max = b[0]->n + b[1]->n;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
235 b[0]->hits = realloc(b[0]->hits, b[0]->max * sizeof(bsw2hit_t));
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
236 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
237 for (i = 0; i < b[1]->n; ++i) {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
238 bsw2hit_t *p = b[0]->hits + b[0]->n + i;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
239 *p = b[1]->hits[i];
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
240 if (is_reverse) {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
241 int x = p->beg;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
242 p->beg = l - p->end;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
243 p->end = l - x;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
244 p->flag |= 0x10;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
245 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
246 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
247 b[0]->n += b[1]->n;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
248 bsw2_destroy(b[1]);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
249 b[1] = 0;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
250 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
251
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
252 static bwtsw2_t *bsw2_aln1_core(const bsw2opt_t *opt, const bntseq_t *bns, uint8_t *pac, const bwt_t *target, int l, uint8_t *seq[2], int is_rev, bsw2global_t *pool)
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
253 {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
254 extern void bsw2_chain_filter(const bsw2opt_t *opt, int len, bwtsw2_t *b[2]);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
255 bwtsw2_t *b[2], **bb[2];
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
256 int k;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
257 for (k = 0; k < 2; ++k) {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
258 bwtl_t *query = bwtl_seq2bwtl(l, seq[k]);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
259 bb[k] = bsw2_core(opt, query, target, pool);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
260 bwtl_destroy(query);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
261 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
262 b[0] = bb[0][1]; b[1] = bb[1][1];
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
263 bsw2_chain_filter(opt, l, b);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
264 for (k = 0; k < 2; ++k) {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
265 bsw2_extend_left(opt, bb[k][1], seq[k], l, pac, bns->l_pac, is_rev, pool->aln_mem);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
266 merge_hits(bb[k], l, 0);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
267 bsw2_resolve_duphits(0, bb[k][0], 0);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
268 bsw2_extend_rght(opt, bb[k][0], seq[k], l, pac, bns->l_pac, is_rev, pool->aln_mem);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
269 b[k] = bb[k][0];
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
270 free(bb[k]);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
271 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
272 merge_hits(b, l, 1);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
273 bsw2_resolve_query_overlaps(b[0], opt->mask_level);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
274 return b[0];
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
275 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
276
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
277 /* set ->flag to records the origin of the hit (to forward bwt or reverse bwt) */
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
278 static void flag_fr(bwtsw2_t *b[2])
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
279 {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
280 int i, j;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
281 for (i = 0; i < b[0]->n; ++i) {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
282 bsw2hit_t *p = b[0]->hits + i;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
283 p->flag |= 0x10000;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
284 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
285 for (i = 0; i < b[1]->n; ++i) {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
286 bsw2hit_t *p = b[1]->hits + i;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
287 p->flag |= 0x20000;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
288 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
289 for (i = 0; i < b[0]->n; ++i) {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
290 bsw2hit_t *p = b[0]->hits + i;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
291 for (j = 0; j < b[1]->n; ++j) {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
292 bsw2hit_t *q = b[1]->hits + i;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
293 if (q->beg == p->beg && q->end == p->end && q->k == p->k && q->len == p->len && q->G == p->G) {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
294 q->flag |= 0x30000; p->flag |= 0x30000;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
295 break;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
296 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
297 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
298 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
299 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
300
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
301 typedef struct {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
302 int l, tid;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
303 char *name, *seq, *qual, *sam;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
304 } bsw2seq1_t;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
305
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
306 typedef struct {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
307 int n, max;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
308 bsw2seq1_t *seq;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
309 } bsw2seq_t;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
310
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
311 #ifdef HAVE_PTHREAD
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
312 static pthread_mutex_t g_dbwtsw_lock = PTHREAD_MUTEX_INITIALIZER;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
313 #endif
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
314
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
315 static int fix_cigar(const char *qname, const bntseq_t *bns, bsw2hit_t *p, int n_cigar, uint32_t *cigar)
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
316 {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
317 // FIXME: this routine does not work if the query bridge three reference sequences
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
318 int32_t coor, refl, lq;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
319 int x, y, i, seqid;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
320 bns_coor_pac2real(bns, p->k, p->len, &seqid);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
321 coor = p->k - bns->anns[seqid].offset;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
322 refl = bns->anns[seqid].len;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
323 x = coor; y = 0;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
324 // test if the alignment goes beyond the boundary
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
325 for (i = 0; i < n_cigar; ++i) {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
326 int op = cigar[i]&0xf, ln = cigar[i]>>4;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
327 if (op == 1 || op == 4 || op == 5) y += ln;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
328 else if (op == 2) x += ln;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
329 else x += ln, y += ln;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
330 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
331 lq = y; // length of the query sequence
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
332 if (x > refl) { // then fix it
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
333 int j, nc, mq[2], nlen[2];
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
334 uint32_t *cn, kk = 0;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
335 nc = mq[0] = mq[1] = nlen[0] = nlen[1] = 0;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
336 cn = calloc(n_cigar + 3, 4);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
337 x = coor; y = 0;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
338 for (i = j = 0; i < n_cigar; ++i) {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
339 int op = cigar[i]&0xf, ln = cigar[i]>>4;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
340 if (op == 4 || op == 5 || op == 1) { // ins or clipping
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
341 y += ln;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
342 cn[j++] = cigar[i];
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
343 } else if (op == 2) { // del
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
344 if (x + ln >= refl && nc == 0) {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
345 cn[j++] = (uint32_t)(lq - y)<<4 | 4;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
346 nc = j;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
347 cn[j++] = (uint32_t)y<<4 | 4;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
348 kk = p->k + (x + ln - refl);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
349 nlen[0] = x - coor;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
350 nlen[1] = p->len - nlen[0] - ln;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
351 } else cn[j++] = cigar[i];
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
352 x += ln;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
353 } else if (op == 0) { // match
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
354 if (x + ln >= refl && nc == 0) {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
355 // FIXME: not consider a special case where a split right between M and I
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
356 cn[j++] = (uint32_t)(refl - x)<<4 | 0; // write M
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
357 cn[j++] = (uint32_t)(lq - y - (refl - x))<<4 | 4; // write S
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
358 nc = j;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
359 mq[0] += refl - x;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
360 cn[j++] = (uint32_t)(y + (refl - x))<<4 | 4;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
361 if (x + ln - refl) cn[j++] = (uint32_t)(x + ln - refl)<<4 | 0;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
362 mq[1] += x + ln - refl;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
363 kk = bns->anns[seqid].offset + refl;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
364 nlen[0] = refl - coor;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
365 nlen[1] = p->len - nlen[0];
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
366 } else {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
367 cn[j++] = cigar[i];
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
368 mq[nc?1:0] += ln;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
369 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
370 x += ln; y += ln;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
371 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
372 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
373 if (mq[0] > mq[1]) { // then take the first alignment
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
374 n_cigar = nc;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
375 memcpy(cigar, cn, 4 * nc);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
376 p->len = nlen[0];
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
377 } else {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
378 p->k = kk; p->len = nlen[1];
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
379 n_cigar = j - nc;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
380 memcpy(cigar, cn + nc, 4 * (j - nc));
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
381 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
382 free(cn);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
383 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
384 return n_cigar;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
385 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
386
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
387 /* generate SAM lines for a sequence in ks with alignment stored in
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
388 * b. ks->name and ks->seq will be freed and set to NULL in the end. */
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
389 static void print_hits(const bntseq_t *bns, const bsw2opt_t *opt, bsw2seq1_t *ks, bwtsw2_t *b)
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
390 {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
391 int i, k;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
392 kstring_t str;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
393 memset(&str, 0, sizeof(kstring_t));
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
394 for (i = 0; i < b->n; ++i) {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
395 bsw2hit_t *p = b->hits + i;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
396 int32_t seqid = -1, coor = -1;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
397 int j, qual, nn = 0;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
398 int beg, end;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
399 if (p->l == 0) {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
400 b->n_cigar[i] = fix_cigar(ks->name, bns, p, b->n_cigar[i], b->cigar[i]);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
401 nn = bns_coor_pac2real(bns, p->k, p->len, &seqid);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
402 coor = p->k - bns->anns[seqid].offset;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
403 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
404 ksprintf(&str, "%s\t%d", ks->name, p->flag&0x10);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
405 ksprintf(&str, "\t%s\t%d", seqid>=0? bns->anns[seqid].name : "*", coor + 1);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
406 if (p->l == 0) {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
407 { // estimate mapping quality
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
408 float c = 1.0;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
409 int subo = p->G2 > opt->t? p->G2 : opt->t;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
410 if (p->flag>>16 == 1 || p->flag>>16 == 2) c *= .5;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
411 if (p->n_seeds < 2) c *= .2;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
412 qual = (int)(c * (p->G - subo) * (250.0 / p->G + 0.03 / opt->a) + .499);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
413 if (qual > 250) qual = 250;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
414 if (p->flag&1) qual = 0;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
415 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
416 ksprintf(&str, "\t%d\t", qual);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
417 for (k = 0; k < b->n_cigar[i]; ++k)
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
418 ksprintf(&str, "%d%c", b->cigar[i][k]>>4, (opt->hard_clip? "MIDNHHP" : "MIDNSHP")[b->cigar[i][k]&0xf]);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
419 } else ksprintf(&str, "\t0\t*");
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
420 ksprintf(&str, "\t*\t0\t0\t");
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
421 beg = 0; end = ks->l;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
422 if (opt->hard_clip) {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
423 if ((b->cigar[i][0]&0xf) == 4) beg += b->cigar[i][0]>>4;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
424 if ((b->cigar[i][b->n_cigar[i]-1]&0xf) == 4) end -= b->cigar[i][b->n_cigar[i]-1]>>4;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
425 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
426 for (j = beg; j < end; ++j) {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
427 if (p->flag&0x10) kputc(nt_comp_table[(int)ks->seq[ks->l - 1 - j]], &str);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
428 else kputc(ks->seq[j], &str);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
429 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
430 if (ks->qual) {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
431 kputc('\t', &str);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
432 for (j = beg; j < end; ++j) {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
433 if (p->flag&0x10) kputc(ks->qual[ks->l - 1 - j], &str);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
434 else kputc(ks->qual[j], &str);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
435 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
436 } else ksprintf(&str, "\t*");
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
437 ksprintf(&str, "\tAS:i:%d\tXS:i:%d\tXF:i:%d\tXE:i:%d\tXN:i:%d", p->G, p->G2, p->flag>>16, p->n_seeds, nn);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
438 if (p->l) ksprintf(&str, "\tXI:i:%d", p->l - p->k + 1);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
439 kputc('\n', &str);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
440 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
441 ks->sam = str.s;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
442 free(ks->seq); ks->seq = 0;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
443 free(ks->qual); ks->qual = 0;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
444 free(ks->name); ks->name = 0;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
445 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
446
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
447 /* Core routine to align reads in _seq. It is separated from
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
448 * process_seqs() to realize multi-threading */
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
449 static void bsw2_aln_core(int tid, bsw2seq_t *_seq, const bsw2opt_t *_opt, const bntseq_t *bns, uint8_t *pac, bwt_t * const target[2])
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
450 {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
451 int x;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
452 bsw2opt_t opt = *_opt;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
453 bsw2global_t *pool = bsw2_global_init();
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
454 for (x = 0; x < _seq->n; ++x) {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
455 bsw2seq1_t *p = _seq->seq + x;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
456 uint8_t *seq[2], *rseq[2];
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
457 int i, l, k;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
458 bwtsw2_t *b[2];
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
459 l = p->l;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
460
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
461 #ifdef HAVE_PTHREAD
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
462 if (_opt->n_threads > 1) {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
463 pthread_mutex_lock(&g_dbwtsw_lock);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
464 if (p->tid < 0) p->tid = tid;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
465 else if (p->tid != tid) {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
466 pthread_mutex_unlock(&g_dbwtsw_lock);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
467 continue;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
468 } // in pinciple else should not happen
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
469 pthread_mutex_unlock(&g_dbwtsw_lock);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
470 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
471 #endif
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
472
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
473 // set opt->t
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
474 opt.t = _opt->t;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
475 if (opt.t < log(l) * opt.coef) opt.t = (int)(log(l) * opt.coef + .499);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
476 if (pool->max_l < l) { // then enlarge working space for aln_extend_core()
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
477 int tmp = ((l + 1) / 2 * opt.a + opt.r) / opt.r + l;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
478 pool->max_l = l;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
479 pool->aln_mem = realloc(pool->aln_mem, (tmp + 2) * 24);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
480 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
481 // set opt->bw
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
482 opt.bw = _opt->bw;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
483 k = (l * opt.a - 2 * opt.q) / (2 * opt.r + opt.a);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
484 i = (l * opt.a - opt.a - opt.t) / opt.r;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
485 if (k > i) k = i;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
486 if (k < 1) k = 1; // I do not know if k==0 causes troubles
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
487 opt.bw = _opt->bw < k? _opt->bw : k;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
488 // set seq[2] and rseq[2]
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
489 seq[0] = calloc(l * 4, 1);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
490 seq[1] = seq[0] + l;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
491 rseq[0] = seq[1] + l; rseq[1] = rseq[0] + l;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
492 // convert sequences to 2-bit representation
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
493 for (i = k = 0; i < l; ++i) {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
494 int c = nst_nt4_table[(int)p->seq[i]];
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
495 if (c >= 4) { c = (int)(drand48() * 4); ++k; } // FIXME: ambiguous bases are not properly handled
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
496 seq[0][i] = c;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
497 seq[1][l-1-i] = 3 - c;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
498 rseq[0][l-1-i] = c;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
499 rseq[1][i] = 3 - c;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
500 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
501 if (l - k < opt.t) {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
502 free(seq[0]); continue;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
503 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
504 // alignment
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
505 b[0] = bsw2_aln1_core(&opt, bns, pac, target[0], l, seq, 0, pool);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
506 for (k = 0; k < b[0]->n; ++k)
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
507 if (b[0]->hits[k].n_seeds < opt.t_seeds) break;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
508 if (k < b[0]->n) {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
509 b[1] = bsw2_aln1_core(&opt, bns, pac, target[1], l, rseq, 1, pool);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
510 for (i = 0; i < b[1]->n; ++i) {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
511 bsw2hit_t *p = b[1]->hits + i;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
512 int x = p->beg;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
513 p->beg = l - p->end;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
514 p->end = l - x;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
515 if (p->l == 0) p->k = bns->l_pac - (p->k + p->len);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
516 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
517 flag_fr(b);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
518 merge_hits(b, l, 0);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
519 bsw2_resolve_duphits(0, b[0], 0);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
520 bsw2_resolve_query_overlaps(b[0], opt.mask_level);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
521 } else b[1] = 0;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
522 // generate CIGAR and print SAM
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
523 gen_cigar(&opt, l, seq, pac, b[0]);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
524 print_hits(bns, &opt, p, b[0]);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
525 // free
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
526 free(seq[0]);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
527 bsw2_destroy(b[0]);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
528 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
529 bsw2_global_destroy(pool);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
530 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
531
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
532 #ifdef HAVE_PTHREAD
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
533 typedef struct {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
534 int tid;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
535 bsw2seq_t *_seq;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
536 const bsw2opt_t *_opt;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
537 const bntseq_t *bns;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
538 uint8_t *pac;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
539 bwt_t *target[2];
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
540 } thread_aux_t;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
541
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
542 /* another interface to bsw2_aln_core() to facilitate pthread_create() */
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
543 static void *worker(void *data)
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
544 {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
545 thread_aux_t *p = (thread_aux_t*)data;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
546 bsw2_aln_core(p->tid, p->_seq, p->_opt, p->bns, p->pac, p->target);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
547 return 0;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
548 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
549 #endif
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
550
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
551 /* process sequences stored in _seq, generate SAM lines for these
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
552 * sequences and reset _seq afterwards. */
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
553 static void process_seqs(bsw2seq_t *_seq, const bsw2opt_t *opt, const bntseq_t *bns, uint8_t *pac, bwt_t * const target[2])
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
554 {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
555 int i;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
556
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
557 #ifdef HAVE_PTHREAD
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
558 if (opt->n_threads <= 1) {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
559 bsw2_aln_core(0, _seq, opt, bns, pac, target);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
560 } else {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
561 pthread_t *tid;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
562 pthread_attr_t attr;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
563 thread_aux_t *data;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
564 int j;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
565 pthread_attr_init(&attr);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
566 pthread_attr_setdetachstate(&attr, PTHREAD_CREATE_JOINABLE);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
567 data = (thread_aux_t*)calloc(opt->n_threads, sizeof(thread_aux_t));
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
568 tid = (pthread_t*)calloc(opt->n_threads, sizeof(pthread_t));
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
569 for (j = 0; j < opt->n_threads; ++j) {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
570 thread_aux_t *p = data + j;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
571 p->tid = j; p->_seq = _seq; p->_opt = opt; p->bns = bns;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
572 p->pac = pac; p->target[0] = target[0]; p->target[1] = target[1];
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
573 pthread_create(&tid[j], &attr, worker, p);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
574 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
575 for (j = 0; j < opt->n_threads; ++j) pthread_join(tid[j], 0);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
576 free(data); free(tid);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
577 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
578 #else
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
579 bsw2_aln_core(0, _seq, opt, bns, pac, target);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
580 #endif
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
581
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
582 // print and reset
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
583 for (i = 0; i < _seq->n; ++i) {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
584 bsw2seq1_t *p = _seq->seq + i;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
585 if (p->sam) printf("%s", p->sam);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
586 free(p->name); free(p->seq); free(p->qual); free(p->sam);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
587 p->tid = -1; p->l = 0;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
588 p->name = p->seq = p->qual = p->sam = 0;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
589 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
590 fflush(stdout);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
591 _seq->n = 0;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
592 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
593
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
594 void bsw2_aln(const bsw2opt_t *opt, const bntseq_t *bns, bwt_t * const target[2], const char *fn)
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
595 {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
596 gzFile fp;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
597 kseq_t *ks;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
598 int l, size = 0;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
599 uint8_t *pac;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
600 bsw2seq_t *_seq;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
601
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
602 pac = calloc(bns->l_pac/4+1, 1);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
603 if (pac == 0) {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
604 fprintf(stderr, "[bsw2_aln] insufficient memory!\n");
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
605 return;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
606 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
607 for (l = 0; l < bns->n_seqs; ++l)
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
608 printf("@SQ\tSN:%s\tLN:%d\n", bns->anns[l].name, bns->anns[l].len);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
609 fread(pac, 1, bns->l_pac/4+1, bns->fp_pac);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
610 fp = xzopen(fn, "r");
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
611 ks = kseq_init(fp);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
612 _seq = calloc(1, sizeof(bsw2seq_t));
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
613 while ((l = kseq_read(ks)) >= 0) {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
614 bsw2seq1_t *p;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
615 if (_seq->n == _seq->max) {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
616 _seq->max = _seq->max? _seq->max<<1 : 1024;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
617 _seq->seq = realloc(_seq->seq, _seq->max * sizeof(bsw2seq1_t));
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
618 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
619 p = &_seq->seq[_seq->n++];
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
620 p->tid = -1;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
621 p->l = l;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
622 p->name = strdup(ks->name.s);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
623 p->seq = strdup(ks->seq.s);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
624 p->qual = ks->qual.l? strdup(ks->qual.s) : 0;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
625 p->sam = 0;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
626 size += l;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
627 if (size > opt->chunk_size) {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
628 fprintf(stderr, "[bsw2_aln] read %d sequences (%d bp)...\n", _seq->n, size);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
629 process_seqs(_seq, opt, bns, pac, target);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
630 size = 0;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
631 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
632 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
633 fprintf(stderr, "[bsw2_aln] read %d sequences (%d bp)...\n", _seq->n, size);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
634 process_seqs(_seq, opt, bns, pac, target);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
635 free(_seq->seq); free(_seq);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
636 kseq_destroy(ks);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
637 gzclose(fp);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
638 free(pac);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
639 }