annotate bwa-0.6.2/bwaseqio.c @ 0:dd1186b11b3b draft

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