Mercurial > repos > ashvark > qiime_1_8_0
comparison bwa-0.6.2/bwaseqio.c @ 2:a294fbfcb1db draft default tip
Uploaded BWA
author | ashvark |
---|---|
date | Fri, 18 Jul 2014 07:55:59 -0400 |
parents | dd1186b11b3b |
children |
comparison
equal
deleted
inserted
replaced
1:a9636dc1e99a | 2:a294fbfcb1db |
---|---|
1 #include <zlib.h> | |
2 #include <ctype.h> | |
3 #include "bwtaln.h" | |
4 #include "utils.h" | |
5 #include "bamlite.h" | |
6 | |
7 #include "kseq.h" | |
8 KSEQ_INIT(gzFile, gzread) | |
9 | |
10 extern unsigned char nst_nt4_table[256]; | |
11 static char bam_nt16_nt4_table[] = { 4, 0, 1, 4, 2, 4, 4, 4, 3, 4, 4, 4, 4, 4, 4, 4 }; | |
12 | |
13 struct __bwa_seqio_t { | |
14 // for BAM input | |
15 int is_bam, which; // 1st bit: read1, 2nd bit: read2, 3rd: SE | |
16 bamFile fp; | |
17 // for fastq input | |
18 kseq_t *ks; | |
19 }; | |
20 | |
21 bwa_seqio_t *bwa_bam_open(const char *fn, int which) | |
22 { | |
23 bwa_seqio_t *bs; | |
24 bam_header_t *h; | |
25 bs = (bwa_seqio_t*)calloc(1, sizeof(bwa_seqio_t)); | |
26 bs->is_bam = 1; | |
27 bs->which = which; | |
28 bs->fp = bam_open(fn, "r"); | |
29 h = bam_header_read(bs->fp); | |
30 bam_header_destroy(h); | |
31 return bs; | |
32 } | |
33 | |
34 bwa_seqio_t *bwa_seq_open(const char *fn) | |
35 { | |
36 gzFile fp; | |
37 bwa_seqio_t *bs; | |
38 bs = (bwa_seqio_t*)calloc(1, sizeof(bwa_seqio_t)); | |
39 fp = xzopen(fn, "r"); | |
40 bs->ks = kseq_init(fp); | |
41 return bs; | |
42 } | |
43 | |
44 void bwa_seq_close(bwa_seqio_t *bs) | |
45 { | |
46 if (bs == 0) return; | |
47 if (bs->is_bam) bam_close(bs->fp); | |
48 else { | |
49 gzclose(bs->ks->f->f); | |
50 kseq_destroy(bs->ks); | |
51 } | |
52 free(bs); | |
53 } | |
54 | |
55 void seq_reverse(int len, ubyte_t *seq, int is_comp) | |
56 { | |
57 int i; | |
58 if (is_comp) { | |
59 for (i = 0; i < len>>1; ++i) { | |
60 char tmp = seq[len-1-i]; | |
61 if (tmp < 4) tmp = 3 - tmp; | |
62 seq[len-1-i] = (seq[i] >= 4)? seq[i] : 3 - seq[i]; | |
63 seq[i] = tmp; | |
64 } | |
65 if (len&1) seq[i] = (seq[i] >= 4)? seq[i] : 3 - seq[i]; | |
66 } else { | |
67 for (i = 0; i < len>>1; ++i) { | |
68 char tmp = seq[len-1-i]; | |
69 seq[len-1-i] = seq[i]; seq[i] = tmp; | |
70 } | |
71 } | |
72 } | |
73 | |
74 int bwa_trim_read(int trim_qual, bwa_seq_t *p) | |
75 { | |
76 int s = 0, l, max = 0, max_l = p->len; | |
77 if (trim_qual < 1 || p->qual == 0) return 0; | |
78 for (l = p->len - 1; l >= BWA_MIN_RDLEN; --l) { | |
79 s += trim_qual - (p->qual[l] - 33); | |
80 if (s < 0) break; | |
81 if (s > max) max = s, max_l = l; | |
82 } | |
83 p->clip_len = p->len = max_l; | |
84 return p->full_len - p->len; | |
85 } | |
86 | |
87 static bwa_seq_t *bwa_read_bam(bwa_seqio_t *bs, int n_needed, int *n, int is_comp, int trim_qual) | |
88 { | |
89 bwa_seq_t *seqs, *p; | |
90 int n_seqs, l, i; | |
91 long n_trimmed = 0, n_tot = 0; | |
92 bam1_t *b; | |
93 | |
94 b = bam_init1(); | |
95 n_seqs = 0; | |
96 seqs = (bwa_seq_t*)calloc(n_needed, sizeof(bwa_seq_t)); | |
97 while (bam_read1(bs->fp, b) >= 0) { | |
98 uint8_t *s, *q; | |
99 int go = 0; | |
100 if ((bs->which & 1) && (b->core.flag & BAM_FREAD1)) go = 1; | |
101 if ((bs->which & 2) && (b->core.flag & BAM_FREAD2)) go = 1; | |
102 if ((bs->which & 4) && !(b->core.flag& BAM_FREAD1) && !(b->core.flag& BAM_FREAD2))go = 1; | |
103 if (go == 0) continue; | |
104 l = b->core.l_qseq; | |
105 p = &seqs[n_seqs++]; | |
106 p->tid = -1; // no assigned to a thread | |
107 p->qual = 0; | |
108 p->full_len = p->clip_len = p->len = l; | |
109 n_tot += p->full_len; | |
110 s = bam1_seq(b); q = bam1_qual(b); | |
111 p->seq = (ubyte_t*)calloc(p->len + 1, 1); | |
112 p->qual = (ubyte_t*)calloc(p->len + 1, 1); | |
113 for (i = 0; i != p->full_len; ++i) { | |
114 p->seq[i] = bam_nt16_nt4_table[(int)bam1_seqi(s, i)]; | |
115 p->qual[i] = q[i] + 33 < 126? q[i] + 33 : 126; | |
116 } | |
117 if (bam1_strand(b)) { // then reverse | |
118 seq_reverse(p->len, p->seq, 1); | |
119 seq_reverse(p->len, p->qual, 0); | |
120 } | |
121 if (trim_qual >= 1) n_trimmed += bwa_trim_read(trim_qual, p); | |
122 p->rseq = (ubyte_t*)calloc(p->full_len, 1); | |
123 memcpy(p->rseq, p->seq, p->len); | |
124 seq_reverse(p->len, p->seq, 0); // *IMPORTANT*: will be reversed back in bwa_refine_gapped() | |
125 seq_reverse(p->len, p->rseq, is_comp); | |
126 p->name = strdup((const char*)bam1_qname(b)); | |
127 if (n_seqs == n_needed) break; | |
128 } | |
129 *n = n_seqs; | |
130 if (n_seqs && trim_qual >= 1) | |
131 fprintf(stderr, "[bwa_read_seq] %.1f%% bases are trimmed.\n", 100.0f * n_trimmed/n_tot); | |
132 if (n_seqs == 0) { | |
133 free(seqs); | |
134 bam_destroy1(b); | |
135 return 0; | |
136 } | |
137 bam_destroy1(b); | |
138 return seqs; | |
139 } | |
140 | |
141 #define BARCODE_LOW_QUAL 13 | |
142 | |
143 bwa_seq_t *bwa_read_seq(bwa_seqio_t *bs, int n_needed, int *n, int mode, int trim_qual) | |
144 { | |
145 bwa_seq_t *seqs, *p; | |
146 kseq_t *seq = bs->ks; | |
147 int n_seqs, l, i, is_comp = mode&BWA_MODE_COMPREAD, is_64 = mode&BWA_MODE_IL13, l_bc = mode>>24; | |
148 long n_trimmed = 0, n_tot = 0; | |
149 | |
150 if (l_bc > BWA_MAX_BCLEN) { | |
151 fprintf(stderr, "[%s] the maximum barcode length is %d.\n", __func__, BWA_MAX_BCLEN); | |
152 return 0; | |
153 } | |
154 if (bs->is_bam) return bwa_read_bam(bs, n_needed, n, is_comp, trim_qual); // l_bc has no effect for BAM input | |
155 n_seqs = 0; | |
156 seqs = (bwa_seq_t*)calloc(n_needed, sizeof(bwa_seq_t)); | |
157 while ((l = kseq_read(seq)) >= 0) { | |
158 if ((mode & BWA_MODE_CFY) && (seq->comment.l != 0)) { | |
159 // skip reads that are marked to be filtered by Casava | |
160 char *s = index(seq->comment.s, ':'); | |
161 if (s && *(++s) == 'Y') { | |
162 continue; | |
163 } | |
164 } | |
165 if (is_64 && seq->qual.l) | |
166 for (i = 0; i < seq->qual.l; ++i) seq->qual.s[i] -= 31; | |
167 if (seq->seq.l <= l_bc) continue; // sequence length equals or smaller than the barcode length | |
168 p = &seqs[n_seqs++]; | |
169 if (l_bc) { // then trim barcode | |
170 for (i = 0; i < l_bc; ++i) | |
171 p->bc[i] = (seq->qual.l && seq->qual.s[i]-33 < BARCODE_LOW_QUAL)? tolower(seq->seq.s[i]) : toupper(seq->seq.s[i]); | |
172 p->bc[i] = 0; | |
173 for (; i < seq->seq.l; ++i) | |
174 seq->seq.s[i - l_bc] = seq->seq.s[i]; | |
175 seq->seq.l -= l_bc; seq->seq.s[seq->seq.l] = 0; | |
176 if (seq->qual.l) { | |
177 for (i = l_bc; i < seq->qual.l; ++i) | |
178 seq->qual.s[i - l_bc] = seq->qual.s[i]; | |
179 seq->qual.l -= l_bc; seq->qual.s[seq->qual.l] = 0; | |
180 } | |
181 l = seq->seq.l; | |
182 } else p->bc[0] = 0; | |
183 p->tid = -1; // no assigned to a thread | |
184 p->qual = 0; | |
185 p->full_len = p->clip_len = p->len = l; | |
186 n_tot += p->full_len; | |
187 p->seq = (ubyte_t*)calloc(p->len, 1); | |
188 for (i = 0; i != p->full_len; ++i) | |
189 p->seq[i] = nst_nt4_table[(int)seq->seq.s[i]]; | |
190 if (seq->qual.l) { // copy quality | |
191 p->qual = (ubyte_t*)strdup((char*)seq->qual.s); | |
192 if (trim_qual >= 1) n_trimmed += bwa_trim_read(trim_qual, p); | |
193 } | |
194 p->rseq = (ubyte_t*)calloc(p->full_len, 1); | |
195 memcpy(p->rseq, p->seq, p->len); | |
196 seq_reverse(p->len, p->seq, 0); // *IMPORTANT*: will be reversed back in bwa_refine_gapped() | |
197 seq_reverse(p->len, p->rseq, is_comp); | |
198 p->name = strdup((const char*)seq->name.s); | |
199 { // trim /[12]$ | |
200 int t = strlen(p->name); | |
201 if (t > 2 && p->name[t-2] == '/' && (p->name[t-1] == '1' || p->name[t-1] == '2')) p->name[t-2] = '\0'; | |
202 } | |
203 if (n_seqs == n_needed) break; | |
204 } | |
205 *n = n_seqs; | |
206 if (n_seqs && trim_qual >= 1) | |
207 fprintf(stderr, "[bwa_read_seq] %.1f%% bases are trimmed.\n", 100.0f * n_trimmed/n_tot); | |
208 if (n_seqs == 0) { | |
209 free(seqs); | |
210 return 0; | |
211 } | |
212 return seqs; | |
213 } | |
214 | |
215 void bwa_free_read_seq(int n_seqs, bwa_seq_t *seqs) | |
216 { | |
217 int i, j; | |
218 for (i = 0; i != n_seqs; ++i) { | |
219 bwa_seq_t *p = seqs + i; | |
220 for (j = 0; j < p->n_multi; ++j) | |
221 if (p->multi[j].cigar) free(p->multi[j].cigar); | |
222 free(p->name); | |
223 free(p->seq); free(p->rseq); free(p->qual); free(p->aln); free(p->md); free(p->multi); | |
224 free(p->cigar); | |
225 } | |
226 free(seqs); | |
227 } |