|
0
|
1 /*
|
|
|
2 Author: petr.danecek@sanger
|
|
|
3 gcc -Wall -Winline -g -O2 -I ~/git/samtools bamcheck.c -o bamcheck -lm -lz -L ~/git/samtools -lbam -lpthread
|
|
|
4
|
|
|
5 Assumptions, approximations and other issues:
|
|
|
6 - GC-depth graph does not split reads, the starting position determines which bin is incremented.
|
|
|
7 There are small overlaps between bins (max readlen-1). However, the bins are big (20k).
|
|
|
8 - coverage distribution ignores softclips and deletions
|
|
|
9 - some stats require sorted BAMs
|
|
|
10 - GC content graph can have an untidy, step-like pattern when BAM contains multiple read lengths.
|
|
|
11 - 'bases mapped' (stats->nbases_mapped) is calculated from read lengths given by BAM (core.l_qseq)
|
|
|
12 - With the -t option, the whole reads are used. Except for the number of mapped bases (cigar)
|
|
|
13 counts, no splicing is done, no indels or soft clips are considered, even small overlap is
|
|
|
14 good enough to include the read in the stats.
|
|
|
15
|
|
|
16 */
|
|
|
17
|
|
|
18 #define BAMCHECK_VERSION "2012-09-04"
|
|
|
19
|
|
|
20 #define _ISOC99_SOURCE
|
|
|
21 #include <stdio.h>
|
|
|
22 #include <stdlib.h>
|
|
|
23 #include <stdarg.h>
|
|
|
24 #include <string.h>
|
|
|
25 #include <math.h>
|
|
|
26 #include <ctype.h>
|
|
|
27 #include <getopt.h>
|
|
|
28 #include <errno.h>
|
|
|
29 #include <assert.h>
|
|
|
30 #include "faidx.h"
|
|
|
31 #include "khash.h"
|
|
|
32 #include "sam.h"
|
|
|
33 #include "sam_header.h"
|
|
|
34 #include "razf.h"
|
|
|
35
|
|
|
36 #define BWA_MIN_RDLEN 35
|
|
|
37 #define IS_PAIRED(bam) ((bam)->core.flag&BAM_FPAIRED && !((bam)->core.flag&BAM_FUNMAP) && !((bam)->core.flag&BAM_FMUNMAP))
|
|
|
38 #define IS_UNMAPPED(bam) ((bam)->core.flag&BAM_FUNMAP)
|
|
|
39 #define IS_REVERSE(bam) ((bam)->core.flag&BAM_FREVERSE)
|
|
|
40 #define IS_MATE_REVERSE(bam) ((bam)->core.flag&BAM_FMREVERSE)
|
|
|
41 #define IS_READ1(bam) ((bam)->core.flag&BAM_FREAD1)
|
|
|
42 #define IS_READ2(bam) ((bam)->core.flag&BAM_FREAD2)
|
|
|
43 #define IS_DUP(bam) ((bam)->core.flag&BAM_FDUP)
|
|
|
44
|
|
|
45 typedef struct
|
|
|
46 {
|
|
|
47 int32_t line_len, line_blen;
|
|
|
48 int64_t len;
|
|
|
49 uint64_t offset;
|
|
|
50 }
|
|
|
51 faidx1_t;
|
|
|
52 KHASH_MAP_INIT_STR(kh_faidx, faidx1_t)
|
|
|
53 KHASH_MAP_INIT_STR(kh_bam_tid, int)
|
|
|
54 KHASH_MAP_INIT_STR(kh_rg, const char *)
|
|
|
55 struct __faidx_t {
|
|
|
56 RAZF *rz;
|
|
|
57 int n, m;
|
|
|
58 char **name;
|
|
|
59 khash_t(kh_faidx) *hash;
|
|
|
60 };
|
|
|
61
|
|
|
62 typedef struct
|
|
|
63 {
|
|
|
64 float gc;
|
|
|
65 uint32_t depth;
|
|
|
66 }
|
|
|
67 gc_depth_t;
|
|
|
68
|
|
|
69 // For coverage distribution, a simple pileup
|
|
|
70 typedef struct
|
|
|
71 {
|
|
|
72 int64_t pos;
|
|
|
73 int size, start;
|
|
|
74 int *buffer;
|
|
|
75 }
|
|
|
76 round_buffer_t;
|
|
|
77
|
|
|
78 typedef struct { uint32_t from, to; } pos_t;
|
|
|
79 typedef struct
|
|
|
80 {
|
|
|
81 int npos,mpos,cpos;
|
|
|
82 pos_t *pos;
|
|
|
83 }
|
|
|
84 regions_t;
|
|
|
85
|
|
|
86 typedef struct
|
|
|
87 {
|
|
|
88 // Parameters
|
|
|
89 int trim_qual; // bwa trim quality
|
|
|
90
|
|
|
91 // Dimensions of the quality histogram holder (quals_1st,quals_2nd), GC content holder (gc_1st,gc_2nd),
|
|
|
92 // insert size histogram holder
|
|
|
93 int nquals; // The number of quality bins
|
|
|
94 int nbases; // The maximum sequence length the allocated array can hold
|
|
|
95 int nisize; // The maximum insert size that the allocated array can hold
|
|
|
96 int ngc; // The size of gc_1st and gc_2nd
|
|
|
97 int nindels; // The maximum indel length for indel distribution
|
|
|
98
|
|
|
99 // Arrays for the histogram data
|
|
|
100 uint64_t *quals_1st, *quals_2nd;
|
|
|
101 uint64_t *gc_1st, *gc_2nd;
|
|
|
102 uint64_t *isize_inward, *isize_outward, *isize_other;
|
|
|
103 uint64_t *acgt_cycles;
|
|
|
104 uint64_t *read_lengths;
|
|
|
105 uint64_t *insertions, *deletions;
|
|
|
106 uint64_t *ins_cycles_1st, *ins_cycles_2nd, *del_cycles_1st, *del_cycles_2nd;
|
|
|
107
|
|
|
108 // The extremes encountered
|
|
|
109 int max_len; // Maximum read length
|
|
|
110 int max_qual; // Maximum quality
|
|
|
111 float isize_main_bulk; // There are always some unrealistically big insert sizes, report only the main part
|
|
|
112 int is_sorted;
|
|
|
113
|
|
|
114 // Summary numbers
|
|
|
115 uint64_t total_len;
|
|
|
116 uint64_t total_len_dup;
|
|
|
117 uint64_t nreads_1st;
|
|
|
118 uint64_t nreads_2nd;
|
|
|
119 uint64_t nreads_filtered;
|
|
|
120 uint64_t nreads_dup;
|
|
|
121 uint64_t nreads_unmapped;
|
|
|
122 uint64_t nreads_unpaired;
|
|
|
123 uint64_t nreads_paired;
|
|
|
124 uint64_t nreads_anomalous;
|
|
|
125 uint64_t nreads_mq0;
|
|
|
126 uint64_t nbases_mapped;
|
|
|
127 uint64_t nbases_mapped_cigar;
|
|
|
128 uint64_t nbases_trimmed; // bwa trimmed bases
|
|
|
129 uint64_t nmismatches;
|
|
|
130 uint64_t nreads_QCfailed, nreads_secondary;
|
|
|
131
|
|
|
132 // GC-depth related data
|
|
|
133 uint32_t ngcd, igcd; // The maximum number of GC depth bins and index of the current bin
|
|
|
134 gc_depth_t *gcd; // The GC-depth bins holder
|
|
|
135 int gcd_bin_size; // The size of GC-depth bin
|
|
|
136 uint32_t gcd_ref_size; // The approximate size of the genome
|
|
|
137 int32_t tid, gcd_pos; // Position of the current bin
|
|
|
138 int32_t pos; // Position of the last read
|
|
|
139
|
|
|
140 // Coverage distribution related data
|
|
|
141 int ncov; // The number of coverage bins
|
|
|
142 uint64_t *cov; // The coverage frequencies
|
|
|
143 int cov_min,cov_max,cov_step; // Minimum, maximum coverage and size of the coverage bins
|
|
|
144 round_buffer_t cov_rbuf; // Pileup round buffer
|
|
|
145
|
|
|
146 // Mismatches by read cycle
|
|
|
147 uint8_t *rseq_buf; // A buffer for reference sequence to check the mismatches against
|
|
|
148 int mrseq_buf; // The size of the buffer
|
|
|
149 int32_t rseq_pos; // The coordinate of the first base in the buffer
|
|
|
150 int32_t nrseq_buf; // The used part of the buffer
|
|
|
151 uint64_t *mpc_buf; // Mismatches per cycle
|
|
|
152
|
|
|
153 // Filters
|
|
|
154 int filter_readlen;
|
|
|
155
|
|
|
156 // Target regions
|
|
|
157 int nregions, reg_from,reg_to;
|
|
|
158 regions_t *regions;
|
|
|
159
|
|
|
160 // Auxiliary data
|
|
|
161 int flag_require, flag_filter;
|
|
|
162 double sum_qual; // For calculating average quality value
|
|
|
163 samfile_t *sam;
|
|
|
164 khash_t(kh_rg) *rg_hash; // Read groups to include, the array is null-terminated
|
|
|
165 faidx_t *fai; // Reference sequence for GC-depth graph
|
|
|
166 int argc; // Command line arguments to be printed on the output
|
|
|
167 char **argv;
|
|
|
168 }
|
|
|
169 stats_t;
|
|
|
170
|
|
|
171 void error(const char *format, ...);
|
|
|
172 void bam_init_header_hash(bam_header_t *header);
|
|
|
173 int is_in_regions(bam1_t *bam_line, stats_t *stats);
|
|
|
174
|
|
|
175
|
|
|
176 // Coverage distribution methods
|
|
|
177 inline int coverage_idx(int min, int max, int n, int step, int depth)
|
|
|
178 {
|
|
|
179 if ( depth < min )
|
|
|
180 return 0;
|
|
|
181
|
|
|
182 if ( depth > max )
|
|
|
183 return n-1;
|
|
|
184
|
|
|
185 return 1 + (depth - min) / step;
|
|
|
186 }
|
|
|
187
|
|
|
188 inline int round_buffer_lidx2ridx(int offset, int size, int64_t refpos, int64_t pos)
|
|
|
189 {
|
|
|
190 return (offset + (pos-refpos) % size) % size;
|
|
|
191 }
|
|
|
192
|
|
|
193 void round_buffer_flush(stats_t *stats, int64_t pos)
|
|
|
194 {
|
|
|
195 int ibuf,idp;
|
|
|
196
|
|
|
197 if ( pos==stats->cov_rbuf.pos )
|
|
|
198 return;
|
|
|
199
|
|
|
200 int64_t new_pos = pos;
|
|
|
201 if ( pos==-1 || pos - stats->cov_rbuf.pos >= stats->cov_rbuf.size )
|
|
|
202 {
|
|
|
203 // Flush the whole buffer, but in sequential order,
|
|
|
204 pos = stats->cov_rbuf.pos + stats->cov_rbuf.size - 1;
|
|
|
205 }
|
|
|
206
|
|
|
207 if ( pos < stats->cov_rbuf.pos )
|
|
|
208 error("Expected coordinates in ascending order, got %ld after %ld\n", pos,stats->cov_rbuf.pos);
|
|
|
209
|
|
|
210 int ifrom = stats->cov_rbuf.start;
|
|
|
211 int ito = round_buffer_lidx2ridx(stats->cov_rbuf.start,stats->cov_rbuf.size,stats->cov_rbuf.pos,pos-1);
|
|
|
212 if ( ifrom>ito )
|
|
|
213 {
|
|
|
214 for (ibuf=ifrom; ibuf<stats->cov_rbuf.size; ibuf++)
|
|
|
215 {
|
|
|
216 if ( !stats->cov_rbuf.buffer[ibuf] )
|
|
|
217 continue;
|
|
|
218 idp = coverage_idx(stats->cov_min,stats->cov_max,stats->ncov,stats->cov_step,stats->cov_rbuf.buffer[ibuf]);
|
|
|
219 stats->cov[idp]++;
|
|
|
220 stats->cov_rbuf.buffer[ibuf] = 0;
|
|
|
221 }
|
|
|
222 ifrom = 0;
|
|
|
223 }
|
|
|
224 for (ibuf=ifrom; ibuf<=ito; ibuf++)
|
|
|
225 {
|
|
|
226 if ( !stats->cov_rbuf.buffer[ibuf] )
|
|
|
227 continue;
|
|
|
228 idp = coverage_idx(stats->cov_min,stats->cov_max,stats->ncov,stats->cov_step,stats->cov_rbuf.buffer[ibuf]);
|
|
|
229 stats->cov[idp]++;
|
|
|
230 stats->cov_rbuf.buffer[ibuf] = 0;
|
|
|
231 }
|
|
|
232 stats->cov_rbuf.start = (new_pos==-1) ? 0 : round_buffer_lidx2ridx(stats->cov_rbuf.start,stats->cov_rbuf.size,stats->cov_rbuf.pos,pos);
|
|
|
233 stats->cov_rbuf.pos = new_pos;
|
|
|
234 }
|
|
|
235
|
|
|
236 void round_buffer_insert_read(round_buffer_t *rbuf, int64_t from, int64_t to)
|
|
|
237 {
|
|
|
238 if ( to-from >= rbuf->size )
|
|
|
239 error("The read length too big (%d), please increase the buffer length (currently %d)\n", to-from+1,rbuf->size);
|
|
|
240 if ( from < rbuf->pos )
|
|
|
241 error("The reads are not sorted (%ld comes after %ld).\n", from,rbuf->pos);
|
|
|
242
|
|
|
243 int ifrom,ito,ibuf;
|
|
|
244 ifrom = round_buffer_lidx2ridx(rbuf->start,rbuf->size,rbuf->pos,from);
|
|
|
245 ito = round_buffer_lidx2ridx(rbuf->start,rbuf->size,rbuf->pos,to);
|
|
|
246 if ( ifrom>ito )
|
|
|
247 {
|
|
|
248 for (ibuf=ifrom; ibuf<rbuf->size; ibuf++)
|
|
|
249 rbuf->buffer[ibuf]++;
|
|
|
250 ifrom = 0;
|
|
|
251 }
|
|
|
252 for (ibuf=ifrom; ibuf<=ito; ibuf++)
|
|
|
253 rbuf->buffer[ibuf]++;
|
|
|
254 }
|
|
|
255
|
|
|
256 // Calculate the number of bases in the read trimmed by BWA
|
|
|
257 int bwa_trim_read(int trim_qual, uint8_t *quals, int len, int reverse)
|
|
|
258 {
|
|
|
259 if ( len<BWA_MIN_RDLEN ) return 0;
|
|
|
260
|
|
|
261 // Although the name implies that the read cannot be trimmed to more than BWA_MIN_RDLEN,
|
|
|
262 // the calculation can in fact trim it to (BWA_MIN_RDLEN-1). (bwa_trim_read in bwa/bwaseqio.c).
|
|
|
263 int max_trimmed = len - BWA_MIN_RDLEN + 1;
|
|
|
264 int l, sum=0, max_sum=0, max_l=0;
|
|
|
265
|
|
|
266 for (l=0; l<max_trimmed; l++)
|
|
|
267 {
|
|
|
268 sum += trim_qual - quals[ reverse ? l : len-1-l ];
|
|
|
269 if ( sum<0 ) break;
|
|
|
270 if ( sum>max_sum )
|
|
|
271 {
|
|
|
272 max_sum = sum;
|
|
|
273 // This is the correct way, but bwa clips from some reason one base less
|
|
|
274 // max_l = l+1;
|
|
|
275 max_l = l;
|
|
|
276 }
|
|
|
277 }
|
|
|
278 return max_l;
|
|
|
279 }
|
|
|
280
|
|
|
281
|
|
|
282 void count_indels(stats_t *stats,bam1_t *bam_line)
|
|
|
283 {
|
|
|
284 int is_fwd = IS_REVERSE(bam_line) ? 0 : 1;
|
|
|
285 int is_1st = IS_READ1(bam_line) ? 1 : 0;
|
|
|
286 int icig;
|
|
|
287 int icycle = 0;
|
|
|
288 int read_len = bam_line->core.l_qseq;
|
|
|
289 for (icig=0; icig<bam_line->core.n_cigar; icig++)
|
|
|
290 {
|
|
|
291 // Conversion from uint32_t to MIDNSHP
|
|
|
292 // 0123456
|
|
|
293 // MIDNSHP
|
|
|
294 int cig = bam1_cigar(bam_line)[icig] & BAM_CIGAR_MASK;
|
|
|
295 int ncig = bam1_cigar(bam_line)[icig] >> BAM_CIGAR_SHIFT;
|
|
|
296
|
|
|
297 if ( cig==1 )
|
|
|
298 {
|
|
|
299 int idx = is_fwd ? icycle : read_len-icycle-ncig;
|
|
|
300 if ( idx<0 )
|
|
|
301 error("FIXME: read_len=%d vs icycle=%d\n", read_len,icycle);
|
|
|
302 if ( idx >= stats->nbases || idx<0 ) error("FIXME: %d vs %d, %s:%d %s\n", idx,stats->nbases, stats->sam->header->target_name[bam_line->core.tid],bam_line->core.pos+1,bam1_qname(bam_line));
|
|
|
303 if ( is_1st )
|
|
|
304 stats->ins_cycles_1st[idx]++;
|
|
|
305 else
|
|
|
306 stats->ins_cycles_2nd[idx]++;
|
|
|
307 icycle += ncig;
|
|
|
308 if ( ncig<=stats->nindels )
|
|
|
309 stats->insertions[ncig-1]++;
|
|
|
310 continue;
|
|
|
311 }
|
|
|
312 if ( cig==2 )
|
|
|
313 {
|
|
|
314 int idx = is_fwd ? icycle-1 : read_len-icycle-1;
|
|
|
315 if ( idx<0 ) continue; // discard meaningless deletions
|
|
|
316 if ( idx >= stats->nbases ) error("FIXME: %d vs %d\n", idx,stats->nbases);
|
|
|
317 if ( is_1st )
|
|
|
318 stats->del_cycles_1st[idx]++;
|
|
|
319 else
|
|
|
320 stats->del_cycles_2nd[idx]++;
|
|
|
321 if ( ncig<=stats->nindels )
|
|
|
322 stats->deletions[ncig-1]++;
|
|
|
323 continue;
|
|
|
324 }
|
|
|
325 if ( cig!=3 && cig!=5 )
|
|
|
326 icycle += ncig;
|
|
|
327 }
|
|
|
328 }
|
|
|
329
|
|
|
330 void count_mismatches_per_cycle(stats_t *stats,bam1_t *bam_line)
|
|
|
331 {
|
|
|
332 int is_fwd = IS_REVERSE(bam_line) ? 0 : 1;
|
|
|
333 int icig,iread=0,icycle=0;
|
|
|
334 int iref = bam_line->core.pos - stats->rseq_pos;
|
|
|
335 int read_len = bam_line->core.l_qseq;
|
|
|
336 uint8_t *read = bam1_seq(bam_line);
|
|
|
337 uint8_t *quals = bam1_qual(bam_line);
|
|
|
338 uint64_t *mpc_buf = stats->mpc_buf;
|
|
|
339 for (icig=0; icig<bam_line->core.n_cigar; icig++)
|
|
|
340 {
|
|
|
341 // Conversion from uint32_t to MIDNSHP
|
|
|
342 // 0123456
|
|
|
343 // MIDNSHP
|
|
|
344 int cig = bam1_cigar(bam_line)[icig] & BAM_CIGAR_MASK;
|
|
|
345 int ncig = bam1_cigar(bam_line)[icig] >> BAM_CIGAR_SHIFT;
|
|
|
346 if ( cig==1 )
|
|
|
347 {
|
|
|
348 iread += ncig;
|
|
|
349 icycle += ncig;
|
|
|
350 continue;
|
|
|
351 }
|
|
|
352 if ( cig==2 )
|
|
|
353 {
|
|
|
354 iref += ncig;
|
|
|
355 continue;
|
|
|
356 }
|
|
|
357 if ( cig==4 )
|
|
|
358 {
|
|
|
359 icycle += ncig;
|
|
|
360 // Soft-clips are present in the sequence, but the position of the read marks a start of non-clipped sequence
|
|
|
361 // iref += ncig;
|
|
|
362 iread += ncig;
|
|
|
363 continue;
|
|
|
364 }
|
|
|
365 if ( cig==5 )
|
|
|
366 {
|
|
|
367 icycle += ncig;
|
|
|
368 continue;
|
|
|
369 }
|
|
|
370 // Ignore H and N CIGARs. The letter are inserted e.g. by TopHat and often require very large
|
|
|
371 // chunk of refseq in memory. Not very frequent and not noticable in the stats.
|
|
|
372 if ( cig==3 || cig==5 ) continue;
|
|
|
373 if ( cig!=0 )
|
|
|
374 error("TODO: cigar %d, %s:%d %s\n", cig,stats->sam->header->target_name[bam_line->core.tid],bam_line->core.pos+1,bam1_qname(bam_line));
|
|
|
375
|
|
|
376 if ( ncig+iref > stats->nrseq_buf )
|
|
|
377 error("FIXME: %d+%d > %d, %s, %s:%d\n",ncig,iref,stats->nrseq_buf, bam1_qname(bam_line),stats->sam->header->target_name[bam_line->core.tid],bam_line->core.pos+1);
|
|
|
378
|
|
|
379 int im;
|
|
|
380 for (im=0; im<ncig; im++)
|
|
|
381 {
|
|
|
382 uint8_t cread = bam1_seqi(read,iread);
|
|
|
383 uint8_t cref = stats->rseq_buf[iref];
|
|
|
384
|
|
|
385 // ---------------15
|
|
|
386 // =ACMGRSVTWYHKDBN
|
|
|
387 if ( cread==15 )
|
|
|
388 {
|
|
|
389 int idx = is_fwd ? icycle : read_len-icycle-1;
|
|
|
390 if ( idx>stats->max_len )
|
|
|
391 error("mpc: %d>%d\n",idx,stats->max_len);
|
|
|
392 idx = idx*stats->nquals;
|
|
|
393 if ( idx>=stats->nquals*stats->nbases )
|
|
|
394 error("FIXME: mpc_buf overflow\n");
|
|
|
395 mpc_buf[idx]++;
|
|
|
396 }
|
|
|
397 else if ( cref && cread && cref!=cread )
|
|
|
398 {
|
|
|
399 uint8_t qual = quals[iread] + 1;
|
|
|
400 if ( qual>=stats->nquals )
|
|
|
401 error("TODO: quality too high %d>=%d (%s %d %s)\n", qual,stats->nquals, stats->sam->header->target_name[bam_line->core.tid],bam_line->core.pos+1,bam1_qname(bam_line));
|
|
|
402
|
|
|
403 int idx = is_fwd ? icycle : read_len-icycle-1;
|
|
|
404 if ( idx>stats->max_len )
|
|
|
405 error("mpc: %d>%d\n",idx,stats->max_len);
|
|
|
406
|
|
|
407 idx = idx*stats->nquals + qual;
|
|
|
408 if ( idx>=stats->nquals*stats->nbases )
|
|
|
409 error("FIXME: mpc_buf overflow\n");
|
|
|
410 mpc_buf[idx]++;
|
|
|
411 }
|
|
|
412
|
|
|
413 iref++;
|
|
|
414 iread++;
|
|
|
415 icycle++;
|
|
|
416 }
|
|
|
417 }
|
|
|
418 }
|
|
|
419
|
|
|
420 void read_ref_seq(stats_t *stats,int32_t tid,int32_t pos)
|
|
|
421 {
|
|
|
422 khash_t(kh_faidx) *h;
|
|
|
423 khiter_t iter;
|
|
|
424 faidx1_t val;
|
|
|
425 char *chr, c;
|
|
|
426 faidx_t *fai = stats->fai;
|
|
|
427
|
|
|
428 h = fai->hash;
|
|
|
429 chr = stats->sam->header->target_name[tid];
|
|
|
430
|
|
|
431 // ID of the sequence name
|
|
|
432 iter = kh_get(kh_faidx, h, chr);
|
|
|
433 if (iter == kh_end(h))
|
|
|
434 error("No such reference sequence [%s]?\n", chr);
|
|
|
435 val = kh_value(h, iter);
|
|
|
436
|
|
|
437 // Check the boundaries
|
|
|
438 if (pos >= val.len)
|
|
|
439 error("Was the bam file mapped with the reference sequence supplied?"
|
|
|
440 " A read mapped beyond the end of the chromosome (%s:%d, chromosome length %d).\n", chr,pos,val.len);
|
|
|
441 int size = stats->mrseq_buf;
|
|
|
442 // The buffer extends beyond the chromosome end. Later the rest will be filled with N's.
|
|
|
443 if (size+pos > val.len) size = val.len-pos;
|
|
|
444
|
|
|
445 // Position the razf reader
|
|
|
446 razf_seek(fai->rz, val.offset + pos / val.line_blen * val.line_len + pos % val.line_blen, SEEK_SET);
|
|
|
447
|
|
|
448 uint8_t *ptr = stats->rseq_buf;
|
|
|
449 int nread = 0;
|
|
|
450 while ( nread<size && razf_read(fai->rz,&c,1) && !fai->rz->z_err )
|
|
|
451 {
|
|
|
452 if ( !isgraph(c) )
|
|
|
453 continue;
|
|
|
454
|
|
|
455 // Conversion between uint8_t coding and ACGT
|
|
|
456 // -12-4---8-------
|
|
|
457 // =ACMGRSVTWYHKDBN
|
|
|
458 if ( c=='A' || c=='a' )
|
|
|
459 *ptr = 1;
|
|
|
460 else if ( c=='C' || c=='c' )
|
|
|
461 *ptr = 2;
|
|
|
462 else if ( c=='G' || c=='g' )
|
|
|
463 *ptr = 4;
|
|
|
464 else if ( c=='T' || c=='t' )
|
|
|
465 *ptr = 8;
|
|
|
466 else
|
|
|
467 *ptr = 0;
|
|
|
468 ptr++;
|
|
|
469 nread++;
|
|
|
470 }
|
|
|
471 if ( nread < stats->mrseq_buf )
|
|
|
472 {
|
|
|
473 memset(ptr,0, stats->mrseq_buf - nread);
|
|
|
474 nread = stats->mrseq_buf;
|
|
|
475 }
|
|
|
476 stats->nrseq_buf = nread;
|
|
|
477 stats->rseq_pos = pos;
|
|
|
478 stats->tid = tid;
|
|
|
479 }
|
|
|
480
|
|
|
481 float fai_gc_content(stats_t *stats, int pos, int len)
|
|
|
482 {
|
|
|
483 uint32_t gc,count,c;
|
|
|
484 int i = pos - stats->rseq_pos, ito = i + len;
|
|
|
485 assert( i>=0 && ito<=stats->nrseq_buf );
|
|
|
486
|
|
|
487 // Count GC content
|
|
|
488 gc = count = 0;
|
|
|
489 for (; i<ito; i++)
|
|
|
490 {
|
|
|
491 c = stats->rseq_buf[i];
|
|
|
492 if ( c==2 || c==4 )
|
|
|
493 {
|
|
|
494 gc++;
|
|
|
495 count++;
|
|
|
496 }
|
|
|
497 else if ( c==1 || c==8 )
|
|
|
498 count++;
|
|
|
499 }
|
|
|
500 return count ? (float)gc/count : 0;
|
|
|
501 }
|
|
|
502
|
|
|
503 void realloc_rseq_buffer(stats_t *stats)
|
|
|
504 {
|
|
|
505 int n = stats->nbases*10;
|
|
|
506 if ( stats->gcd_bin_size > n ) n = stats->gcd_bin_size;
|
|
|
507 if ( stats->mrseq_buf<n )
|
|
|
508 {
|
|
|
509 stats->rseq_buf = realloc(stats->rseq_buf,sizeof(uint8_t)*n);
|
|
|
510 stats->mrseq_buf = n;
|
|
|
511 }
|
|
|
512 }
|
|
|
513
|
|
|
514 void realloc_gcd_buffer(stats_t *stats, int seq_len)
|
|
|
515 {
|
|
|
516 if ( seq_len >= stats->gcd_bin_size )
|
|
|
517 error("The --GC-depth bin size (%d) is set too low for the read length %d\n", stats->gcd_bin_size, seq_len);
|
|
|
518
|
|
|
519 int n = 1 + stats->gcd_ref_size / (stats->gcd_bin_size - seq_len);
|
|
|
520 if ( n <= stats->igcd )
|
|
|
521 error("The --GC-depth bin size is too small or reference genome too big; please decrease the bin size or increase the reference length\n");
|
|
|
522
|
|
|
523 if ( n > stats->ngcd )
|
|
|
524 {
|
|
|
525 stats->gcd = realloc(stats->gcd, n*sizeof(gc_depth_t));
|
|
|
526 if ( !stats->gcd )
|
|
|
527 error("Could not realloc GCD buffer, too many chromosomes or the genome too long?? [%u %u]\n", stats->ngcd,n);
|
|
|
528 memset(&(stats->gcd[stats->ngcd]),0,(n-stats->ngcd)*sizeof(gc_depth_t));
|
|
|
529 stats->ngcd = n;
|
|
|
530 }
|
|
|
531
|
|
|
532 realloc_rseq_buffer(stats);
|
|
|
533 }
|
|
|
534
|
|
|
535 void realloc_buffers(stats_t *stats, int seq_len)
|
|
|
536 {
|
|
|
537 int n = 2*(1 + seq_len - stats->nbases) + stats->nbases;
|
|
|
538
|
|
|
539 stats->quals_1st = realloc(stats->quals_1st, n*stats->nquals*sizeof(uint64_t));
|
|
|
540 if ( !stats->quals_1st )
|
|
|
541 error("Could not realloc buffers, the sequence too long: %d (%ld)\n", seq_len,n*stats->nquals*sizeof(uint64_t));
|
|
|
542 memset(stats->quals_1st + stats->nbases*stats->nquals, 0, (n-stats->nbases)*stats->nquals*sizeof(uint64_t));
|
|
|
543
|
|
|
544 stats->quals_2nd = realloc(stats->quals_2nd, n*stats->nquals*sizeof(uint64_t));
|
|
|
545 if ( !stats->quals_2nd )
|
|
|
546 error("Could not realloc buffers, the sequence too long: %d (2x%ld)\n", seq_len,n*stats->nquals*sizeof(uint64_t));
|
|
|
547 memset(stats->quals_2nd + stats->nbases*stats->nquals, 0, (n-stats->nbases)*stats->nquals*sizeof(uint64_t));
|
|
|
548
|
|
|
549 if ( stats->mpc_buf )
|
|
|
550 {
|
|
|
551 stats->mpc_buf = realloc(stats->mpc_buf, n*stats->nquals*sizeof(uint64_t));
|
|
|
552 if ( !stats->mpc_buf )
|
|
|
553 error("Could not realloc buffers, the sequence too long: %d (%ld)\n", seq_len,n*stats->nquals*sizeof(uint64_t));
|
|
|
554 memset(stats->mpc_buf + stats->nbases*stats->nquals, 0, (n-stats->nbases)*stats->nquals*sizeof(uint64_t));
|
|
|
555 }
|
|
|
556
|
|
|
557 stats->acgt_cycles = realloc(stats->acgt_cycles, n*4*sizeof(uint64_t));
|
|
|
558 if ( !stats->acgt_cycles )
|
|
|
559 error("Could not realloc buffers, the sequence too long: %d (%ld)\n", seq_len,n*4*sizeof(uint64_t));
|
|
|
560 memset(stats->acgt_cycles + stats->nbases*4, 0, (n-stats->nbases)*4*sizeof(uint64_t));
|
|
|
561
|
|
|
562 stats->read_lengths = realloc(stats->read_lengths, n*sizeof(uint64_t));
|
|
|
563 if ( !stats->read_lengths )
|
|
|
564 error("Could not realloc buffers, the sequence too long: %d (%ld)\n", seq_len,n*sizeof(uint64_t));
|
|
|
565 memset(stats->read_lengths + stats->nbases, 0, (n-stats->nbases)*sizeof(uint64_t));
|
|
|
566
|
|
|
567 stats->insertions = realloc(stats->insertions, n*sizeof(uint64_t));
|
|
|
568 if ( !stats->insertions )
|
|
|
569 error("Could not realloc buffers, the sequence too long: %d (%ld)\n", seq_len,n*sizeof(uint64_t));
|
|
|
570 memset(stats->insertions + stats->nbases, 0, (n-stats->nbases)*sizeof(uint64_t));
|
|
|
571
|
|
|
572 stats->deletions = realloc(stats->deletions, n*sizeof(uint64_t));
|
|
|
573 if ( !stats->deletions )
|
|
|
574 error("Could not realloc buffers, the sequence too long: %d (%ld)\n", seq_len,n*sizeof(uint64_t));
|
|
|
575 memset(stats->deletions + stats->nbases, 0, (n-stats->nbases)*sizeof(uint64_t));
|
|
|
576
|
|
|
577 stats->ins_cycles_1st = realloc(stats->ins_cycles_1st, (n+1)*sizeof(uint64_t));
|
|
|
578 if ( !stats->ins_cycles_1st )
|
|
|
579 error("Could not realloc buffers, the sequence too long: %d (%ld)\n", seq_len,(n+1)*sizeof(uint64_t));
|
|
|
580 memset(stats->ins_cycles_1st + stats->nbases + 1, 0, (n-stats->nbases)*sizeof(uint64_t));
|
|
|
581
|
|
|
582 stats->ins_cycles_2nd = realloc(stats->ins_cycles_2nd, (n+1)*sizeof(uint64_t));
|
|
|
583 if ( !stats->ins_cycles_2nd )
|
|
|
584 error("Could not realloc buffers, the sequence too long: %d (%ld)\n", seq_len,(n+1)*sizeof(uint64_t));
|
|
|
585 memset(stats->ins_cycles_2nd + stats->nbases + 1, 0, (n-stats->nbases)*sizeof(uint64_t));
|
|
|
586
|
|
|
587 stats->del_cycles_1st = realloc(stats->del_cycles_1st, (n+1)*sizeof(uint64_t));
|
|
|
588 if ( !stats->del_cycles_1st )
|
|
|
589 error("Could not realloc buffers, the sequence too long: %d (%ld)\n", seq_len,(n+1)*sizeof(uint64_t));
|
|
|
590 memset(stats->del_cycles_1st + stats->nbases + 1, 0, (n-stats->nbases)*sizeof(uint64_t));
|
|
|
591
|
|
|
592 stats->del_cycles_2nd = realloc(stats->del_cycles_2nd, (n+1)*sizeof(uint64_t));
|
|
|
593 if ( !stats->del_cycles_2nd )
|
|
|
594 error("Could not realloc buffers, the sequence too long: %d (%ld)\n", seq_len,(n+1)*sizeof(uint64_t));
|
|
|
595 memset(stats->del_cycles_2nd + stats->nbases + 1, 0, (n-stats->nbases)*sizeof(uint64_t));
|
|
|
596
|
|
|
597 stats->nbases = n;
|
|
|
598
|
|
|
599 // Realloc the coverage distribution buffer
|
|
|
600 int *rbuffer = calloc(sizeof(int),seq_len*5);
|
|
|
601 n = stats->cov_rbuf.size-stats->cov_rbuf.start;
|
|
|
602 memcpy(rbuffer,stats->cov_rbuf.buffer+stats->cov_rbuf.start,n);
|
|
|
603 if ( stats->cov_rbuf.start>1 )
|
|
|
604 memcpy(rbuffer+n,stats->cov_rbuf.buffer,stats->cov_rbuf.start);
|
|
|
605 stats->cov_rbuf.start = 0;
|
|
|
606 free(stats->cov_rbuf.buffer);
|
|
|
607 stats->cov_rbuf.buffer = rbuffer;
|
|
|
608 stats->cov_rbuf.size = seq_len*5;
|
|
|
609
|
|
|
610 realloc_rseq_buffer(stats);
|
|
|
611 }
|
|
|
612
|
|
|
613 void collect_stats(bam1_t *bam_line, stats_t *stats)
|
|
|
614 {
|
|
|
615 if ( stats->rg_hash )
|
|
|
616 {
|
|
|
617 const uint8_t *rg = bam_aux_get(bam_line, "RG");
|
|
|
618 if ( !rg ) return;
|
|
|
619 khiter_t k = kh_get(kh_rg, stats->rg_hash, (const char*)(rg + 1));
|
|
|
620 if ( k == kh_end(stats->rg_hash) ) return;
|
|
|
621 }
|
|
|
622 if ( stats->flag_require && (bam_line->core.flag & stats->flag_require)!=stats->flag_require )
|
|
|
623 {
|
|
|
624 stats->nreads_filtered++;
|
|
|
625 return;
|
|
|
626 }
|
|
|
627 if ( stats->flag_filter && (bam_line->core.flag & stats->flag_filter) )
|
|
|
628 {
|
|
|
629 stats->nreads_filtered++;
|
|
|
630 return;
|
|
|
631 }
|
|
|
632 if ( !is_in_regions(bam_line,stats) )
|
|
|
633 return;
|
|
|
634 if ( stats->filter_readlen!=-1 && bam_line->core.l_qseq!=stats->filter_readlen )
|
|
|
635 return;
|
|
|
636
|
|
|
637 if ( bam_line->core.flag & BAM_FQCFAIL ) stats->nreads_QCfailed++;
|
|
|
638 if ( bam_line->core.flag & BAM_FSECONDARY ) stats->nreads_secondary++;
|
|
|
639
|
|
|
640 int seq_len = bam_line->core.l_qseq;
|
|
|
641 if ( !seq_len ) return;
|
|
|
642
|
|
|
643 if ( seq_len >= stats->nbases )
|
|
|
644 realloc_buffers(stats,seq_len);
|
|
|
645 if ( stats->max_len<seq_len )
|
|
|
646 stats->max_len = seq_len;
|
|
|
647
|
|
|
648 stats->read_lengths[seq_len]++;
|
|
|
649
|
|
|
650 // Count GC and ACGT per cycle
|
|
|
651 uint8_t base, *seq = bam1_seq(bam_line);
|
|
|
652 int gc_count = 0;
|
|
|
653 int i;
|
|
|
654 int reverse = IS_REVERSE(bam_line);
|
|
|
655 for (i=0; i<seq_len; i++)
|
|
|
656 {
|
|
|
657 // Conversion from uint8_t coding to ACGT
|
|
|
658 // -12-4---8-------
|
|
|
659 // =ACMGRSVTWYHKDBN
|
|
|
660 // 01 2 3
|
|
|
661 base = bam1_seqi(seq,i);
|
|
|
662 base /= 2;
|
|
|
663 if ( base==1 || base==2 ) gc_count++;
|
|
|
664 else if ( base>2 ) base=3;
|
|
|
665 if ( 4*(reverse ? seq_len-i-1 : i) + base >= stats->nbases*4 )
|
|
|
666 error("FIXME: acgt_cycles\n");
|
|
|
667 stats->acgt_cycles[ 4*(reverse ? seq_len-i-1 : i) + base ]++;
|
|
|
668 }
|
|
|
669 int gc_idx_min = gc_count*(stats->ngc-1)/seq_len;
|
|
|
670 int gc_idx_max = (gc_count+1)*(stats->ngc-1)/seq_len;
|
|
|
671 if ( gc_idx_max >= stats->ngc ) gc_idx_max = stats->ngc - 1;
|
|
|
672
|
|
|
673 // Determine which array (1st or 2nd read) will these stats go to,
|
|
|
674 // trim low quality bases from end the same way BWA does,
|
|
|
675 // fill GC histogram
|
|
|
676 uint64_t *quals;
|
|
|
677 uint8_t *bam_quals = bam1_qual(bam_line);
|
|
|
678 if ( bam_line->core.flag&BAM_FREAD2 )
|
|
|
679 {
|
|
|
680 quals = stats->quals_2nd;
|
|
|
681 stats->nreads_2nd++;
|
|
|
682 for (i=gc_idx_min; i<gc_idx_max; i++)
|
|
|
683 stats->gc_2nd[i]++;
|
|
|
684 }
|
|
|
685 else
|
|
|
686 {
|
|
|
687 quals = stats->quals_1st;
|
|
|
688 stats->nreads_1st++;
|
|
|
689 for (i=gc_idx_min; i<gc_idx_max; i++)
|
|
|
690 stats->gc_1st[i]++;
|
|
|
691 }
|
|
|
692 if ( stats->trim_qual>0 )
|
|
|
693 stats->nbases_trimmed += bwa_trim_read(stats->trim_qual, bam_quals, seq_len, reverse);
|
|
|
694
|
|
|
695 // Quality histogram and average quality
|
|
|
696 for (i=0; i<seq_len; i++)
|
|
|
697 {
|
|
|
698 uint8_t qual = bam_quals[ reverse ? seq_len-i-1 : i];
|
|
|
699 if ( qual>=stats->nquals )
|
|
|
700 error("TODO: quality too high %d>=%d (%s %d %s)\n", qual,stats->nquals,stats->sam->header->target_name[bam_line->core.tid],bam_line->core.pos+1,bam1_qname(bam_line));
|
|
|
701 if ( qual>stats->max_qual )
|
|
|
702 stats->max_qual = qual;
|
|
|
703
|
|
|
704 quals[ i*stats->nquals+qual ]++;
|
|
|
705 stats->sum_qual += qual;
|
|
|
706 }
|
|
|
707
|
|
|
708 // Look at the flags and increment appropriate counters (mapped, paired, etc)
|
|
|
709 if ( IS_UNMAPPED(bam_line) )
|
|
|
710 stats->nreads_unmapped++;
|
|
|
711 else
|
|
|
712 {
|
|
|
713 if ( !bam_line->core.qual )
|
|
|
714 stats->nreads_mq0++;
|
|
|
715
|
|
|
716 count_indels(stats,bam_line);
|
|
|
717
|
|
|
718 if ( !IS_PAIRED(bam_line) )
|
|
|
719 stats->nreads_unpaired++;
|
|
|
720 else
|
|
|
721 {
|
|
|
722 stats->nreads_paired++;
|
|
|
723
|
|
|
724 if ( bam_line->core.tid!=bam_line->core.mtid )
|
|
|
725 stats->nreads_anomalous++;
|
|
|
726
|
|
|
727 // The insert size is tricky, because for long inserts the libraries are
|
|
|
728 // prepared differently and the pairs point in other direction. BWA does
|
|
|
729 // not set the paired flag for them. Similar thing is true also for 454
|
|
|
730 // reads. Mates mapped to different chromosomes have isize==0.
|
|
|
731 int32_t isize = bam_line->core.isize;
|
|
|
732 if ( isize<0 ) isize = -isize;
|
|
|
733 if ( isize >= stats->nisize )
|
|
|
734 isize = stats->nisize-1;
|
|
|
735 if ( isize>0 || bam_line->core.tid==bam_line->core.mtid )
|
|
|
736 {
|
|
|
737 int pos_fst = bam_line->core.mpos - bam_line->core.pos;
|
|
|
738 int is_fst = IS_READ1(bam_line) ? 1 : -1;
|
|
|
739 int is_fwd = IS_REVERSE(bam_line) ? -1 : 1;
|
|
|
740 int is_mfwd = IS_MATE_REVERSE(bam_line) ? -1 : 1;
|
|
|
741
|
|
|
742 if ( is_fwd*is_mfwd>0 )
|
|
|
743 stats->isize_other[isize]++;
|
|
|
744 else if ( is_fst*pos_fst>0 )
|
|
|
745 {
|
|
|
746 if ( is_fst*is_fwd>0 )
|
|
|
747 stats->isize_inward[isize]++;
|
|
|
748 else
|
|
|
749 stats->isize_outward[isize]++;
|
|
|
750 }
|
|
|
751 else if ( is_fst*pos_fst<0 )
|
|
|
752 {
|
|
|
753 if ( is_fst*is_fwd>0 )
|
|
|
754 stats->isize_outward[isize]++;
|
|
|
755 else
|
|
|
756 stats->isize_inward[isize]++;
|
|
|
757 }
|
|
|
758 }
|
|
|
759 }
|
|
|
760
|
|
|
761 // Number of mismatches
|
|
|
762 uint8_t *nm = bam_aux_get(bam_line,"NM");
|
|
|
763 if (nm)
|
|
|
764 stats->nmismatches += bam_aux2i(nm);
|
|
|
765
|
|
|
766 // Number of mapped bases from cigar
|
|
|
767 // Conversion from uint32_t to MIDNSHP
|
|
|
768 // 012-4--
|
|
|
769 // MIDNSHP
|
|
|
770 if ( bam_line->core.n_cigar == 0)
|
|
|
771 error("FIXME: mapped read with no cigar?\n");
|
|
|
772 int readlen=seq_len;
|
|
|
773 if ( stats->regions )
|
|
|
774 {
|
|
|
775 // Count only on-target bases
|
|
|
776 int iref = bam_line->core.pos + 1;
|
|
|
777 for (i=0; i<bam_line->core.n_cigar; i++)
|
|
|
778 {
|
|
|
779 int cig = bam1_cigar(bam_line)[i]&BAM_CIGAR_MASK;
|
|
|
780 int ncig = bam1_cigar(bam_line)[i]>>BAM_CIGAR_SHIFT;
|
|
|
781 if ( cig==2 ) readlen += ncig;
|
|
|
782 else if ( cig==0 )
|
|
|
783 {
|
|
|
784 if ( iref < stats->reg_from ) ncig -= stats->reg_from-iref;
|
|
|
785 else if ( iref+ncig-1 > stats->reg_to ) ncig -= iref+ncig-1 - stats->reg_to;
|
|
|
786 if ( ncig<0 ) ncig = 0;
|
|
|
787 stats->nbases_mapped_cigar += ncig;
|
|
|
788 iref += bam1_cigar(bam_line)[i]>>BAM_CIGAR_SHIFT;
|
|
|
789 }
|
|
|
790 else if ( cig==1 )
|
|
|
791 {
|
|
|
792 iref += ncig;
|
|
|
793 if ( iref>=stats->reg_from && iref<=stats->reg_to )
|
|
|
794 stats->nbases_mapped_cigar += ncig;
|
|
|
795 }
|
|
|
796 }
|
|
|
797 }
|
|
|
798 else
|
|
|
799 {
|
|
|
800 // Count the whole read
|
|
|
801 for (i=0; i<bam_line->core.n_cigar; i++)
|
|
|
802 {
|
|
|
803 if ( (bam1_cigar(bam_line)[i]&BAM_CIGAR_MASK)==0 || (bam1_cigar(bam_line)[i]&BAM_CIGAR_MASK)==1 )
|
|
|
804 stats->nbases_mapped_cigar += bam1_cigar(bam_line)[i]>>BAM_CIGAR_SHIFT;
|
|
|
805 if ( (bam1_cigar(bam_line)[i]&BAM_CIGAR_MASK)==2 )
|
|
|
806 readlen += bam1_cigar(bam_line)[i]>>BAM_CIGAR_SHIFT;
|
|
|
807 }
|
|
|
808 }
|
|
|
809 stats->nbases_mapped += seq_len;
|
|
|
810
|
|
|
811 if ( stats->tid==bam_line->core.tid && bam_line->core.pos<stats->pos )
|
|
|
812 stats->is_sorted = 0;
|
|
|
813 stats->pos = bam_line->core.pos;
|
|
|
814
|
|
|
815 if ( stats->is_sorted )
|
|
|
816 {
|
|
|
817 if ( stats->tid==-1 || stats->tid!=bam_line->core.tid )
|
|
|
818 round_buffer_flush(stats,-1);
|
|
|
819
|
|
|
820 // Mismatches per cycle and GC-depth graph. For simplicity, reads overlapping GCD bins
|
|
|
821 // are not splitted which results in up to seq_len-1 overlaps. The default bin size is
|
|
|
822 // 20kbp, so the effect is negligible.
|
|
|
823 if ( stats->fai )
|
|
|
824 {
|
|
|
825 int inc_ref = 0, inc_gcd = 0;
|
|
|
826 // First pass or new chromosome
|
|
|
827 if ( stats->rseq_pos==-1 || stats->tid != bam_line->core.tid ) { inc_ref=1; inc_gcd=1; }
|
|
|
828 // Read goes beyond the end of the rseq buffer
|
|
|
829 else if ( stats->rseq_pos+stats->nrseq_buf < bam_line->core.pos+readlen ) { inc_ref=1; inc_gcd=1; }
|
|
|
830 // Read overlaps the next gcd bin
|
|
|
831 else if ( stats->gcd_pos+stats->gcd_bin_size < bam_line->core.pos+readlen )
|
|
|
832 {
|
|
|
833 inc_gcd = 1;
|
|
|
834 if ( stats->rseq_pos+stats->nrseq_buf < bam_line->core.pos+stats->gcd_bin_size ) inc_ref = 1;
|
|
|
835 }
|
|
|
836 if ( inc_gcd )
|
|
|
837 {
|
|
|
838 stats->igcd++;
|
|
|
839 if ( stats->igcd >= stats->ngcd )
|
|
|
840 realloc_gcd_buffer(stats, readlen);
|
|
|
841 if ( inc_ref )
|
|
|
842 read_ref_seq(stats,bam_line->core.tid,bam_line->core.pos);
|
|
|
843 stats->gcd_pos = bam_line->core.pos;
|
|
|
844 stats->gcd[ stats->igcd ].gc = fai_gc_content(stats, stats->gcd_pos, stats->gcd_bin_size);
|
|
|
845 }
|
|
|
846
|
|
|
847 count_mismatches_per_cycle(stats,bam_line);
|
|
|
848 }
|
|
|
849 // No reference and first pass, new chromosome or sequence going beyond the end of the gcd bin
|
|
|
850 else if ( stats->gcd_pos==-1 || stats->tid != bam_line->core.tid || bam_line->core.pos - stats->gcd_pos > stats->gcd_bin_size )
|
|
|
851 {
|
|
|
852 // First pass or a new chromosome
|
|
|
853 stats->tid = bam_line->core.tid;
|
|
|
854 stats->gcd_pos = bam_line->core.pos;
|
|
|
855 stats->igcd++;
|
|
|
856 if ( stats->igcd >= stats->ngcd )
|
|
|
857 realloc_gcd_buffer(stats, readlen);
|
|
|
858 }
|
|
|
859 stats->gcd[ stats->igcd ].depth++;
|
|
|
860 // When no reference sequence is given, approximate the GC from the read (much shorter window, but otherwise OK)
|
|
|
861 if ( !stats->fai )
|
|
|
862 stats->gcd[ stats->igcd ].gc += (float) gc_count / seq_len;
|
|
|
863
|
|
|
864 // Coverage distribution graph
|
|
|
865 round_buffer_flush(stats,bam_line->core.pos);
|
|
|
866 round_buffer_insert_read(&(stats->cov_rbuf),bam_line->core.pos,bam_line->core.pos+seq_len-1);
|
|
|
867 }
|
|
|
868 }
|
|
|
869
|
|
|
870 stats->total_len += seq_len;
|
|
|
871 if ( IS_DUP(bam_line) )
|
|
|
872 {
|
|
|
873 stats->total_len_dup += seq_len;
|
|
|
874 stats->nreads_dup++;
|
|
|
875 }
|
|
|
876 }
|
|
|
877
|
|
|
878 // Sort by GC and depth
|
|
|
879 #define GCD_t(x) ((gc_depth_t *)x)
|
|
|
880 static int gcd_cmp(const void *a, const void *b)
|
|
|
881 {
|
|
|
882 if ( GCD_t(a)->gc < GCD_t(b)->gc ) return -1;
|
|
|
883 if ( GCD_t(a)->gc > GCD_t(b)->gc ) return 1;
|
|
|
884 if ( GCD_t(a)->depth < GCD_t(b)->depth ) return -1;
|
|
|
885 if ( GCD_t(a)->depth > GCD_t(b)->depth ) return 1;
|
|
|
886 return 0;
|
|
|
887 }
|
|
|
888 #undef GCD_t
|
|
|
889
|
|
|
890 float gcd_percentile(gc_depth_t *gcd, int N, int p)
|
|
|
891 {
|
|
|
892 float n,d;
|
|
|
893 int k;
|
|
|
894
|
|
|
895 n = p*(N+1)/100;
|
|
|
896 k = n;
|
|
|
897 if ( k<=0 )
|
|
|
898 return gcd[0].depth;
|
|
|
899 if ( k>=N )
|
|
|
900 return gcd[N-1].depth;
|
|
|
901
|
|
|
902 d = n - k;
|
|
|
903 return gcd[k-1].depth + d*(gcd[k].depth - gcd[k-1].depth);
|
|
|
904 }
|
|
|
905
|
|
|
906 void output_stats(stats_t *stats)
|
|
|
907 {
|
|
|
908 // Calculate average insert size and standard deviation (from the main bulk data only)
|
|
|
909 int isize, ibulk=0;
|
|
|
910 uint64_t nisize=0, nisize_inward=0, nisize_outward=0, nisize_other=0;
|
|
|
911 for (isize=0; isize<stats->nisize; isize++)
|
|
|
912 {
|
|
|
913 // Each pair was counted twice
|
|
|
914 stats->isize_inward[isize] *= 0.5;
|
|
|
915 stats->isize_outward[isize] *= 0.5;
|
|
|
916 stats->isize_other[isize] *= 0.5;
|
|
|
917
|
|
|
918 nisize_inward += stats->isize_inward[isize];
|
|
|
919 nisize_outward += stats->isize_outward[isize];
|
|
|
920 nisize_other += stats->isize_other[isize];
|
|
|
921 nisize += stats->isize_inward[isize] + stats->isize_outward[isize] + stats->isize_other[isize];
|
|
|
922 }
|
|
|
923
|
|
|
924 double bulk=0, avg_isize=0, sd_isize=0;
|
|
|
925 for (isize=0; isize<stats->nisize; isize++)
|
|
|
926 {
|
|
|
927 bulk += stats->isize_inward[isize] + stats->isize_outward[isize] + stats->isize_other[isize];
|
|
|
928 avg_isize += isize * (stats->isize_inward[isize] + stats->isize_outward[isize] + stats->isize_other[isize]);
|
|
|
929
|
|
|
930 if ( bulk/nisize > stats->isize_main_bulk )
|
|
|
931 {
|
|
|
932 ibulk = isize+1;
|
|
|
933 nisize = bulk;
|
|
|
934 break;
|
|
|
935 }
|
|
|
936 }
|
|
|
937 avg_isize /= nisize ? nisize : 1;
|
|
|
938 for (isize=1; isize<ibulk; isize++)
|
|
|
939 sd_isize += (stats->isize_inward[isize] + stats->isize_outward[isize] + stats->isize_other[isize]) * (isize-avg_isize)*(isize-avg_isize) / nisize;
|
|
|
940 sd_isize = sqrt(sd_isize);
|
|
|
941
|
|
|
942
|
|
|
943 printf("# This file was produced by bamcheck (%s)\n",BAMCHECK_VERSION);
|
|
|
944 printf("# The command line was: %s",stats->argv[0]);
|
|
|
945 int i;
|
|
|
946 for (i=1; i<stats->argc; i++)
|
|
|
947 printf(" %s",stats->argv[i]);
|
|
|
948 printf("\n");
|
|
|
949 printf("# Summary Numbers. Use `grep ^SN | cut -f 2-` to extract this part.\n");
|
|
|
950 printf("SN\traw total sequences:\t%ld\n", (long)(stats->nreads_filtered+stats->nreads_1st+stats->nreads_2nd));
|
|
|
951 printf("SN\tfiltered sequences:\t%ld\n", (long)stats->nreads_filtered);
|
|
|
952 printf("SN\tsequences:\t%ld\n", (long)(stats->nreads_1st+stats->nreads_2nd));
|
|
|
953 printf("SN\tis paired:\t%d\n", stats->nreads_1st&&stats->nreads_2nd ? 1 : 0);
|
|
|
954 printf("SN\tis sorted:\t%d\n", stats->is_sorted ? 1 : 0);
|
|
|
955 printf("SN\t1st fragments:\t%ld\n", (long)stats->nreads_1st);
|
|
|
956 printf("SN\tlast fragments:\t%ld\n", (long)stats->nreads_2nd);
|
|
|
957 printf("SN\treads mapped:\t%ld\n", (long)(stats->nreads_paired+stats->nreads_unpaired));
|
|
|
958 printf("SN\treads unmapped:\t%ld\n", (long)stats->nreads_unmapped);
|
|
|
959 printf("SN\treads unpaired:\t%ld\n", (long)stats->nreads_unpaired);
|
|
|
960 printf("SN\treads paired:\t%ld\n", (long)stats->nreads_paired);
|
|
|
961 printf("SN\treads duplicated:\t%ld\n", (long)stats->nreads_dup);
|
|
|
962 printf("SN\treads MQ0:\t%ld\n", (long)stats->nreads_mq0);
|
|
|
963 printf("SN\treads QC failed:\t%ld\n", (long)stats->nreads_QCfailed);
|
|
|
964 printf("SN\tnon-primary alignments:\t%ld\n", (long)stats->nreads_secondary);
|
|
|
965 printf("SN\ttotal length:\t%ld\n", (long)stats->total_len);
|
|
|
966 printf("SN\tbases mapped:\t%ld\n", (long)stats->nbases_mapped);
|
|
|
967 printf("SN\tbases mapped (cigar):\t%ld\n", (long)stats->nbases_mapped_cigar);
|
|
|
968 printf("SN\tbases trimmed:\t%ld\n", (long)stats->nbases_trimmed);
|
|
|
969 printf("SN\tbases duplicated:\t%ld\n", (long)stats->total_len_dup);
|
|
|
970 printf("SN\tmismatches:\t%ld\n", (long)stats->nmismatches);
|
|
|
971 printf("SN\terror rate:\t%e\n", (float)stats->nmismatches/stats->nbases_mapped_cigar);
|
|
|
972 float avg_read_length = (stats->nreads_1st+stats->nreads_2nd)?stats->total_len/(stats->nreads_1st+stats->nreads_2nd):0;
|
|
|
973 printf("SN\taverage length:\t%.0f\n", avg_read_length);
|
|
|
974 printf("SN\tmaximum length:\t%d\n", stats->max_len);
|
|
|
975 printf("SN\taverage quality:\t%.1f\n", stats->total_len?stats->sum_qual/stats->total_len:0);
|
|
|
976 printf("SN\tinsert size average:\t%.1f\n", avg_isize);
|
|
|
977 printf("SN\tinsert size standard deviation:\t%.1f\n", sd_isize);
|
|
|
978 printf("SN\tinward oriented pairs:\t%ld\n", (long)nisize_inward);
|
|
|
979 printf("SN\toutward oriented pairs:\t%ld\n", (long)nisize_outward);
|
|
|
980 printf("SN\tpairs with other orientation:\t%ld\n", (long)nisize_other);
|
|
|
981 printf("SN\tpairs on different chromosomes:\t%ld\n", (long)stats->nreads_anomalous/2);
|
|
|
982
|
|
|
983 int ibase,iqual;
|
|
|
984 if ( stats->max_len<stats->nbases ) stats->max_len++;
|
|
|
985 if ( stats->max_qual+1<stats->nquals ) stats->max_qual++;
|
|
|
986 printf("# First Fragment Qualitites. Use `grep ^FFQ | cut -f 2-` to extract this part.\n");
|
|
|
987 printf("# Columns correspond to qualities and rows to cycles. First column is the cycle number.\n");
|
|
|
988 for (ibase=0; ibase<stats->max_len; ibase++)
|
|
|
989 {
|
|
|
990 printf("FFQ\t%d",ibase+1);
|
|
|
991 for (iqual=0; iqual<=stats->max_qual; iqual++)
|
|
|
992 {
|
|
|
993 printf("\t%ld", (long)stats->quals_1st[ibase*stats->nquals+iqual]);
|
|
|
994 }
|
|
|
995 printf("\n");
|
|
|
996 }
|
|
|
997 printf("# Last Fragment Qualitites. Use `grep ^LFQ | cut -f 2-` to extract this part.\n");
|
|
|
998 printf("# Columns correspond to qualities and rows to cycles. First column is the cycle number.\n");
|
|
|
999 for (ibase=0; ibase<stats->max_len; ibase++)
|
|
|
1000 {
|
|
|
1001 printf("LFQ\t%d",ibase+1);
|
|
|
1002 for (iqual=0; iqual<=stats->max_qual; iqual++)
|
|
|
1003 {
|
|
|
1004 printf("\t%ld", (long)stats->quals_2nd[ibase*stats->nquals+iqual]);
|
|
|
1005 }
|
|
|
1006 printf("\n");
|
|
|
1007 }
|
|
|
1008 if ( stats->mpc_buf )
|
|
|
1009 {
|
|
|
1010 printf("# Mismatches per cycle and quality. Use `grep ^MPC | cut -f 2-` to extract this part.\n");
|
|
|
1011 printf("# Columns correspond to qualities, rows to cycles. First column is the cycle number, second\n");
|
|
|
1012 printf("# is the number of N's and the rest is the number of mismatches\n");
|
|
|
1013 for (ibase=0; ibase<stats->max_len; ibase++)
|
|
|
1014 {
|
|
|
1015 printf("MPC\t%d",ibase+1);
|
|
|
1016 for (iqual=0; iqual<=stats->max_qual; iqual++)
|
|
|
1017 {
|
|
|
1018 printf("\t%ld", (long)stats->mpc_buf[ibase*stats->nquals+iqual]);
|
|
|
1019 }
|
|
|
1020 printf("\n");
|
|
|
1021 }
|
|
|
1022 }
|
|
|
1023 printf("# GC Content of first fragments. Use `grep ^GCF | cut -f 2-` to extract this part.\n");
|
|
|
1024 int ibase_prev = 0;
|
|
|
1025 for (ibase=0; ibase<stats->ngc; ibase++)
|
|
|
1026 {
|
|
|
1027 if ( stats->gc_1st[ibase]==stats->gc_1st[ibase_prev] ) continue;
|
|
|
1028 printf("GCF\t%.2f\t%ld\n", (ibase+ibase_prev)*0.5*100./(stats->ngc-1), (long)stats->gc_1st[ibase_prev]);
|
|
|
1029 ibase_prev = ibase;
|
|
|
1030 }
|
|
|
1031 printf("# GC Content of last fragments. Use `grep ^GCL | cut -f 2-` to extract this part.\n");
|
|
|
1032 ibase_prev = 0;
|
|
|
1033 for (ibase=0; ibase<stats->ngc; ibase++)
|
|
|
1034 {
|
|
|
1035 if ( stats->gc_2nd[ibase]==stats->gc_2nd[ibase_prev] ) continue;
|
|
|
1036 printf("GCL\t%.2f\t%ld\n", (ibase+ibase_prev)*0.5*100./(stats->ngc-1), (long)stats->gc_2nd[ibase_prev]);
|
|
|
1037 ibase_prev = ibase;
|
|
|
1038 }
|
|
|
1039 printf("# ACGT content per cycle. Use `grep ^GCC | cut -f 2-` to extract this part. The columns are: cycle, and A,C,G,T counts [%%]\n");
|
|
|
1040 for (ibase=0; ibase<stats->max_len; ibase++)
|
|
|
1041 {
|
|
|
1042 uint64_t *ptr = &(stats->acgt_cycles[ibase*4]);
|
|
|
1043 uint64_t sum = ptr[0]+ptr[1]+ptr[2]+ptr[3];
|
|
|
1044 if ( ! sum ) continue;
|
|
|
1045 printf("GCC\t%d\t%.2f\t%.2f\t%.2f\t%.2f\n", ibase,100.*ptr[0]/sum,100.*ptr[1]/sum,100.*ptr[2]/sum,100.*ptr[3]/sum);
|
|
|
1046 }
|
|
|
1047 printf("# Insert sizes. Use `grep ^IS | cut -f 2-` to extract this part. The columns are: pairs total, inward oriented pairs, outward oriented pairs, other pairs\n");
|
|
|
1048 for (isize=0; isize<ibulk; isize++)
|
|
|
1049 printf("IS\t%d\t%ld\t%ld\t%ld\t%ld\n", isize, (long)(stats->isize_inward[isize]+stats->isize_outward[isize]+stats->isize_other[isize]),
|
|
|
1050 (long)stats->isize_inward[isize], (long)stats->isize_outward[isize], (long)stats->isize_other[isize]);
|
|
|
1051
|
|
|
1052 printf("# Read lengths. Use `grep ^RL | cut -f 2-` to extract this part. The columns are: read length, count\n");
|
|
|
1053 int ilen;
|
|
|
1054 for (ilen=0; ilen<stats->max_len; ilen++)
|
|
|
1055 {
|
|
|
1056 if ( stats->read_lengths[ilen]>0 )
|
|
|
1057 printf("RL\t%d\t%ld\n", ilen, (long)stats->read_lengths[ilen]);
|
|
|
1058 }
|
|
|
1059
|
|
|
1060 printf("# Indel distribution. Use `grep ^ID | cut -f 2-` to extract this part. The columns are: length, number of insertions, number of deletions\n");
|
|
|
1061 for (ilen=0; ilen<stats->nindels; ilen++)
|
|
|
1062 {
|
|
|
1063 if ( stats->insertions[ilen]>0 || stats->deletions[ilen]>0 )
|
|
|
1064 printf("ID\t%d\t%ld\t%ld\n", ilen+1, (long)stats->insertions[ilen], (long)stats->deletions[ilen]);
|
|
|
1065 }
|
|
|
1066
|
|
|
1067 printf("# Indels per cycle. Use `grep ^IC | cut -f 2-` to extract this part. The columns are: cycle, number of insertions (fwd), .. (rev) , number of deletions (fwd), .. (rev)\n");
|
|
|
1068 for (ilen=0; ilen<=stats->nbases; ilen++)
|
|
|
1069 {
|
|
|
1070 // For deletions we print the index of the cycle before the deleted base (1-based) and for insertions
|
|
|
1071 // the index of the cycle of the first inserted base (also 1-based)
|
|
|
1072 if ( stats->ins_cycles_1st[ilen]>0 || stats->ins_cycles_2nd[ilen]>0 || stats->del_cycles_1st[ilen]>0 || stats->del_cycles_2nd[ilen]>0 )
|
|
|
1073 printf("IC\t%d\t%ld\t%ld\t%ld\t%ld\n", ilen+1, (long)stats->ins_cycles_1st[ilen], (long)stats->ins_cycles_2nd[ilen], (long)stats->del_cycles_1st[ilen], (long)stats->del_cycles_2nd[ilen]);
|
|
|
1074 }
|
|
|
1075
|
|
|
1076 printf("# Coverage distribution. Use `grep ^COV | cut -f 2-` to extract this part.\n");
|
|
|
1077 if ( stats->cov[0] )
|
|
|
1078 printf("COV\t[<%d]\t%d\t%ld\n",stats->cov_min,stats->cov_min-1, (long)stats->cov[0]);
|
|
|
1079 int icov;
|
|
|
1080 for (icov=1; icov<stats->ncov-1; icov++)
|
|
|
1081 if ( stats->cov[icov] )
|
|
|
1082 printf("COV\t[%d-%d]\t%d\t%ld\n",stats->cov_min + (icov-1)*stats->cov_step, stats->cov_min + icov*stats->cov_step-1,stats->cov_min + icov*stats->cov_step-1, (long)stats->cov[icov]);
|
|
|
1083 if ( stats->cov[stats->ncov-1] )
|
|
|
1084 printf("COV\t[%d<]\t%d\t%ld\n",stats->cov_min + (stats->ncov-2)*stats->cov_step-1,stats->cov_min + (stats->ncov-2)*stats->cov_step-1, (long)stats->cov[stats->ncov-1]);
|
|
|
1085
|
|
|
1086 // Calculate average GC content, then sort by GC and depth
|
|
|
1087 printf("# GC-depth. Use `grep ^GCD | cut -f 2-` to extract this part. The columns are: GC%%, unique sequence percentiles, 10th, 25th, 50th, 75th and 90th depth percentile\n");
|
|
|
1088 uint32_t igcd;
|
|
|
1089 for (igcd=0; igcd<stats->igcd; igcd++)
|
|
|
1090 {
|
|
|
1091 if ( stats->fai )
|
|
|
1092 stats->gcd[igcd].gc = round(100. * stats->gcd[igcd].gc);
|
|
|
1093 else
|
|
|
1094 if ( stats->gcd[igcd].depth )
|
|
|
1095 stats->gcd[igcd].gc = round(100. * stats->gcd[igcd].gc / stats->gcd[igcd].depth);
|
|
|
1096 }
|
|
|
1097 qsort(stats->gcd, stats->igcd+1, sizeof(gc_depth_t), gcd_cmp);
|
|
|
1098 igcd = 0;
|
|
|
1099 while ( igcd < stats->igcd )
|
|
|
1100 {
|
|
|
1101 // Calculate percentiles (10,25,50,75,90th) for the current GC content and print
|
|
|
1102 uint32_t nbins=0, itmp=igcd;
|
|
|
1103 float gc = stats->gcd[igcd].gc;
|
|
|
1104 while ( itmp<stats->igcd && fabs(stats->gcd[itmp].gc-gc)<0.1 )
|
|
|
1105 {
|
|
|
1106 nbins++;
|
|
|
1107 itmp++;
|
|
|
1108 }
|
|
|
1109 printf("GCD\t%.1f\t%.3f\t%.3f\t%.3f\t%.3f\t%.3f\t%.3f\n", gc, (igcd+nbins+1)*100./(stats->igcd+1),
|
|
|
1110 gcd_percentile(&(stats->gcd[igcd]),nbins,10) *avg_read_length/stats->gcd_bin_size,
|
|
|
1111 gcd_percentile(&(stats->gcd[igcd]),nbins,25) *avg_read_length/stats->gcd_bin_size,
|
|
|
1112 gcd_percentile(&(stats->gcd[igcd]),nbins,50) *avg_read_length/stats->gcd_bin_size,
|
|
|
1113 gcd_percentile(&(stats->gcd[igcd]),nbins,75) *avg_read_length/stats->gcd_bin_size,
|
|
|
1114 gcd_percentile(&(stats->gcd[igcd]),nbins,90) *avg_read_length/stats->gcd_bin_size
|
|
|
1115 );
|
|
|
1116 igcd += nbins;
|
|
|
1117 }
|
|
|
1118 }
|
|
|
1119
|
|
|
1120 size_t mygetline(char **line, size_t *n, FILE *fp)
|
|
|
1121 {
|
|
|
1122 if (line == NULL || n == NULL || fp == NULL)
|
|
|
1123 {
|
|
|
1124 errno = EINVAL;
|
|
|
1125 return -1;
|
|
|
1126 }
|
|
|
1127 if (*n==0 || !*line)
|
|
|
1128 {
|
|
|
1129 *line = NULL;
|
|
|
1130 *n = 0;
|
|
|
1131 }
|
|
|
1132
|
|
|
1133 size_t nread=0;
|
|
|
1134 int c;
|
|
|
1135 while ((c=getc(fp))!= EOF && c!='\n')
|
|
|
1136 {
|
|
|
1137 if ( ++nread>=*n )
|
|
|
1138 {
|
|
|
1139 *n += 255;
|
|
|
1140 *line = realloc(*line, sizeof(char)*(*n));
|
|
|
1141 }
|
|
|
1142 (*line)[nread-1] = c;
|
|
|
1143 }
|
|
|
1144 if ( nread>=*n )
|
|
|
1145 {
|
|
|
1146 *n += 255;
|
|
|
1147 *line = realloc(*line, sizeof(char)*(*n));
|
|
|
1148 }
|
|
|
1149 (*line)[nread] = 0;
|
|
|
1150 return nread>0 ? nread : -1;
|
|
|
1151
|
|
|
1152 }
|
|
|
1153
|
|
|
1154 void init_regions(stats_t *stats, char *file)
|
|
|
1155 {
|
|
|
1156 khiter_t iter;
|
|
|
1157 khash_t(kh_bam_tid) *header_hash;
|
|
|
1158
|
|
|
1159 bam_init_header_hash(stats->sam->header);
|
|
|
1160 header_hash = (khash_t(kh_bam_tid)*)stats->sam->header->hash;
|
|
|
1161
|
|
|
1162 FILE *fp = fopen(file,"r");
|
|
|
1163 if ( !fp ) error("%s: %s\n",file,strerror(errno));
|
|
|
1164
|
|
|
1165 char *line = NULL;
|
|
|
1166 size_t len = 0;
|
|
|
1167 ssize_t nread;
|
|
|
1168 int warned = 0;
|
|
|
1169 int prev_tid=-1, prev_pos=-1;
|
|
|
1170 while ((nread = mygetline(&line, &len, fp)) != -1)
|
|
|
1171 {
|
|
|
1172 if ( line[0] == '#' ) continue;
|
|
|
1173
|
|
|
1174 int i = 0;
|
|
|
1175 while ( i<nread && !isspace(line[i]) ) i++;
|
|
|
1176 if ( i>=nread ) error("Could not parse the file: %s [%s]\n", file,line);
|
|
|
1177 line[i] = 0;
|
|
|
1178
|
|
|
1179 iter = kh_get(kh_bam_tid, header_hash, line);
|
|
|
1180 int tid = kh_val(header_hash, iter);
|
|
|
1181 if ( iter == kh_end(header_hash) )
|
|
|
1182 {
|
|
|
1183 if ( !warned )
|
|
|
1184 fprintf(stderr,"Warning: Some sequences not present in the BAM, e.g. \"%s\". This message is printed only once.\n", line);
|
|
|
1185 warned = 1;
|
|
|
1186 continue;
|
|
|
1187 }
|
|
|
1188
|
|
|
1189 if ( tid >= stats->nregions )
|
|
|
1190 {
|
|
|
1191 stats->regions = realloc(stats->regions,sizeof(regions_t)*(stats->nregions+100));
|
|
|
1192 int j;
|
|
|
1193 for (j=stats->nregions; j<stats->nregions+100; j++)
|
|
|
1194 {
|
|
|
1195 stats->regions[j].npos = stats->regions[j].mpos = stats->regions[j].cpos = 0;
|
|
|
1196 stats->regions[j].pos = NULL;
|
|
|
1197 }
|
|
|
1198 stats->nregions += 100;
|
|
|
1199 }
|
|
|
1200 int npos = stats->regions[tid].npos;
|
|
|
1201 if ( npos >= stats->regions[tid].mpos )
|
|
|
1202 {
|
|
|
1203 stats->regions[tid].mpos += 1000;
|
|
|
1204 stats->regions[tid].pos = realloc(stats->regions[tid].pos,sizeof(pos_t)*stats->regions[tid].mpos);
|
|
|
1205 }
|
|
|
1206
|
|
|
1207 if ( (sscanf(line+i+1,"%d %d",&stats->regions[tid].pos[npos].from,&stats->regions[tid].pos[npos].to))!=2 ) error("Could not parse the region [%s]\n");
|
|
|
1208 if ( prev_tid==-1 || prev_tid!=tid )
|
|
|
1209 {
|
|
|
1210 prev_tid = tid;
|
|
|
1211 prev_pos = stats->regions[tid].pos[npos].from;
|
|
|
1212 }
|
|
|
1213 if ( prev_pos>stats->regions[tid].pos[npos].from )
|
|
|
1214 error("The positions are not in chromosomal order (%s:%d comes after %d)\n", line,stats->regions[tid].pos[npos].from,prev_pos);
|
|
|
1215 stats->regions[tid].npos++;
|
|
|
1216 }
|
|
|
1217 if (line) free(line);
|
|
|
1218 if ( !stats->regions ) error("Unable to map the -t sequences to the BAM sequences.\n");
|
|
|
1219 fclose(fp);
|
|
|
1220 }
|
|
|
1221
|
|
|
1222 void destroy_regions(stats_t *stats)
|
|
|
1223 {
|
|
|
1224 int i;
|
|
|
1225 for (i=0; i<stats->nregions; i++)
|
|
|
1226 {
|
|
|
1227 if ( !stats->regions[i].mpos ) continue;
|
|
|
1228 free(stats->regions[i].pos);
|
|
|
1229 }
|
|
|
1230 if ( stats->regions ) free(stats->regions);
|
|
|
1231 }
|
|
|
1232
|
|
|
1233 static int fetch_read(const bam1_t *bam_line, void *data)
|
|
|
1234 {
|
|
|
1235 collect_stats((bam1_t*)bam_line,(stats_t*)data);
|
|
|
1236 return 1;
|
|
|
1237 }
|
|
|
1238
|
|
|
1239 void reset_regions(stats_t *stats)
|
|
|
1240 {
|
|
|
1241 int i;
|
|
|
1242 for (i=0; i<stats->nregions; i++)
|
|
|
1243 stats->regions[i].cpos = 0;
|
|
|
1244 }
|
|
|
1245
|
|
|
1246 int is_in_regions(bam1_t *bam_line, stats_t *stats)
|
|
|
1247 {
|
|
|
1248 if ( !stats->regions ) return 1;
|
|
|
1249
|
|
|
1250 if ( bam_line->core.tid >= stats->nregions || bam_line->core.tid<0 ) return 0;
|
|
|
1251 if ( !stats->is_sorted ) error("The BAM must be sorted in order for -t to work.\n");
|
|
|
1252
|
|
|
1253 regions_t *reg = &stats->regions[bam_line->core.tid];
|
|
|
1254 if ( reg->cpos==reg->npos ) return 0; // done for this chr
|
|
|
1255
|
|
|
1256 // Find a matching interval or skip this read. No splicing of reads is done, no indels or soft clips considered,
|
|
|
1257 // even small overlap is enough to include the read in the stats.
|
|
|
1258 int i = reg->cpos;
|
|
|
1259 while ( i<reg->npos && reg->pos[i].to<=bam_line->core.pos ) i++;
|
|
|
1260 if ( i>=reg->npos ) { reg->cpos = reg->npos; return 0; }
|
|
|
1261 if ( bam_line->core.pos + bam_line->core.l_qseq + 1 < reg->pos[i].from ) return 0;
|
|
|
1262 reg->cpos = i;
|
|
|
1263 stats->reg_from = reg->pos[i].from;
|
|
|
1264 stats->reg_to = reg->pos[i].to;
|
|
|
1265
|
|
|
1266 return 1;
|
|
|
1267 }
|
|
|
1268
|
|
|
1269 void init_group_id(stats_t *stats, char *id)
|
|
|
1270 {
|
|
|
1271 if ( !stats->sam->header->dict )
|
|
|
1272 stats->sam->header->dict = sam_header_parse2(stats->sam->header->text);
|
|
|
1273 void *iter = stats->sam->header->dict;
|
|
|
1274 const char *key, *val;
|
|
|
1275 int n = 0;
|
|
|
1276 stats->rg_hash = kh_init(kh_rg);
|
|
|
1277 while ( (iter = sam_header2key_val(iter, "RG","ID","SM", &key, &val)) )
|
|
|
1278 {
|
|
|
1279 if ( !strcmp(id,key) || (val && !strcmp(id,val)) )
|
|
|
1280 {
|
|
|
1281 khiter_t k = kh_get(kh_rg, stats->rg_hash, key);
|
|
|
1282 if ( k != kh_end(stats->rg_hash) )
|
|
|
1283 fprintf(stderr, "[init_group_id] The group ID not unique: \"%s\"\n", key);
|
|
|
1284 int ret;
|
|
|
1285 k = kh_put(kh_rg, stats->rg_hash, key, &ret);
|
|
|
1286 kh_value(stats->rg_hash, k) = val;
|
|
|
1287 n++;
|
|
|
1288 }
|
|
|
1289 }
|
|
|
1290 if ( !n )
|
|
|
1291 error("The sample or read group \"%s\" not present.\n", id);
|
|
|
1292 }
|
|
|
1293
|
|
|
1294
|
|
|
1295 void error(const char *format, ...)
|
|
|
1296 {
|
|
|
1297 if ( !format )
|
|
|
1298 {
|
|
|
1299 printf("Version: %s\n", BAMCHECK_VERSION);
|
|
|
1300 printf("About: The program collects statistics from BAM files. The output can be visualized using plot-bamcheck.\n");
|
|
|
1301 printf("Usage: bamcheck [OPTIONS] file.bam\n");
|
|
|
1302 printf(" bamcheck [OPTIONS] file.bam chr:from-to\n");
|
|
|
1303 printf("Options:\n");
|
|
|
1304 printf(" -c, --coverage <int>,<int>,<int> Coverage distribution min,max,step [1,1000,1]\n");
|
|
|
1305 printf(" -d, --remove-dups Exlude from statistics reads marked as duplicates\n");
|
|
|
1306 printf(" -f, --required-flag <int> Required flag, 0 for unset [0]\n");
|
|
|
1307 printf(" -F, --filtering-flag <int> Filtering flag, 0 for unset [0]\n");
|
|
|
1308 printf(" --GC-depth <float,float> Bin size for GC-depth graph and the maximum reference length [2e4,4.2e9]\n");
|
|
|
1309 printf(" -h, --help This help message\n");
|
|
|
1310 printf(" -i, --insert-size <int> Maximum insert size [8000]\n");
|
|
|
1311 printf(" -I, --id <string> Include only listed read group or sample name\n");
|
|
|
1312 printf(" -l, --read-length <int> Include in the statistics only reads with the given read length []\n");
|
|
|
1313 printf(" -m, --most-inserts <float> Report only the main part of inserts [0.99]\n");
|
|
|
1314 printf(" -q, --trim-quality <int> The BWA trimming parameter [0]\n");
|
|
|
1315 printf(" -r, --ref-seq <file> Reference sequence (required for GC-depth calculation).\n");
|
|
|
1316 printf(" -t, --target-regions <file> Do stats in these regions only. Tab-delimited file chr,from,to, 1-based, inclusive.\n");
|
|
|
1317 printf(" -s, --sam Input is SAM\n");
|
|
|
1318 printf("\n");
|
|
|
1319 }
|
|
|
1320 else
|
|
|
1321 {
|
|
|
1322 va_list ap;
|
|
|
1323 va_start(ap, format);
|
|
|
1324 vfprintf(stderr, format, ap);
|
|
|
1325 va_end(ap);
|
|
|
1326 }
|
|
|
1327 exit(-1);
|
|
|
1328 }
|
|
|
1329
|
|
|
1330 int main(int argc, char *argv[])
|
|
|
1331 {
|
|
|
1332 char *targets = NULL;
|
|
|
1333 char *bam_fname = NULL;
|
|
|
1334 char *group_id = NULL;
|
|
|
1335 samfile_t *sam = NULL;
|
|
|
1336 char in_mode[5];
|
|
|
1337
|
|
|
1338 stats_t *stats = calloc(1,sizeof(stats_t));
|
|
|
1339 stats->ngc = 200;
|
|
|
1340 stats->nquals = 256;
|
|
|
1341 stats->nbases = 300;
|
|
|
1342 stats->nisize = 8000;
|
|
|
1343 stats->max_len = 30;
|
|
|
1344 stats->max_qual = 40;
|
|
|
1345 stats->isize_main_bulk = 0.99; // There are always outliers at the far end
|
|
|
1346 stats->gcd_bin_size = 20e3;
|
|
|
1347 stats->gcd_ref_size = 4.2e9;
|
|
|
1348 stats->rseq_pos = -1;
|
|
|
1349 stats->tid = stats->gcd_pos = -1;
|
|
|
1350 stats->igcd = 0;
|
|
|
1351 stats->is_sorted = 1;
|
|
|
1352 stats->cov_min = 1;
|
|
|
1353 stats->cov_max = 1000;
|
|
|
1354 stats->cov_step = 1;
|
|
|
1355 stats->argc = argc;
|
|
|
1356 stats->argv = argv;
|
|
|
1357 stats->filter_readlen = -1;
|
|
|
1358 stats->nindels = stats->nbases;
|
|
|
1359
|
|
|
1360 strcpy(in_mode, "rb");
|
|
|
1361
|
|
|
1362 static struct option loptions[] =
|
|
|
1363 {
|
|
|
1364 {"help",0,0,'h'},
|
|
|
1365 {"remove-dups",0,0,'d'},
|
|
|
1366 {"sam",0,0,'s'},
|
|
|
1367 {"ref-seq",1,0,'r'},
|
|
|
1368 {"coverage",1,0,'c'},
|
|
|
1369 {"read-length",1,0,'l'},
|
|
|
1370 {"insert-size",1,0,'i'},
|
|
|
1371 {"most-inserts",1,0,'m'},
|
|
|
1372 {"trim-quality",1,0,'q'},
|
|
|
1373 {"target-regions",0,0,'t'},
|
|
|
1374 {"required-flag",1,0,'f'},
|
|
|
1375 {"filtering-flag",0,0,'F'},
|
|
|
1376 {"id",1,0,'I'},
|
|
|
1377 {"GC-depth",1,0,1},
|
|
|
1378 {0,0,0,0}
|
|
|
1379 };
|
|
|
1380 int opt;
|
|
|
1381 while ( (opt=getopt_long(argc,argv,"?hdsr:c:l:i:t:m:q:f:F:I:1:",loptions,NULL))>0 )
|
|
|
1382 {
|
|
|
1383 switch (opt)
|
|
|
1384 {
|
|
|
1385 case 'f': stats->flag_require=strtol(optarg,0,0); break;
|
|
|
1386 case 'F': stats->flag_filter=strtol(optarg,0,0); break;
|
|
|
1387 case 'd': stats->flag_filter|=BAM_FDUP; break;
|
|
|
1388 case 's': strcpy(in_mode, "r"); break;
|
|
|
1389 case 'r': stats->fai = fai_load(optarg);
|
|
|
1390 if (stats->fai==0)
|
|
|
1391 error("Could not load faidx: %s\n", optarg);
|
|
|
1392 break;
|
|
|
1393 case 1 : {
|
|
|
1394 float flen,fbin;
|
|
|
1395 if ( sscanf(optarg,"%f,%f",&fbin,&flen)!= 2 )
|
|
|
1396 error("Unable to parse --GC-depth %s\n", optarg);
|
|
|
1397 stats->gcd_bin_size = fbin;
|
|
|
1398 stats->gcd_ref_size = flen;
|
|
|
1399 }
|
|
|
1400 break;
|
|
|
1401 case 'c': if ( sscanf(optarg,"%d,%d,%d",&stats->cov_min,&stats->cov_max,&stats->cov_step)!= 3 )
|
|
|
1402 error("Unable to parse -c %s\n", optarg);
|
|
|
1403 break;
|
|
|
1404 case 'l': stats->filter_readlen = atoi(optarg); break;
|
|
|
1405 case 'i': stats->nisize = atoi(optarg); break;
|
|
|
1406 case 'm': stats->isize_main_bulk = atof(optarg); break;
|
|
|
1407 case 'q': stats->trim_qual = atoi(optarg); break;
|
|
|
1408 case 't': targets = optarg; break;
|
|
|
1409 case 'I': group_id = optarg; break;
|
|
|
1410 case '?':
|
|
|
1411 case 'h': error(NULL);
|
|
|
1412 default: error("Unknown argument: %s\n", optarg);
|
|
|
1413 }
|
|
|
1414 }
|
|
|
1415 if ( optind<argc )
|
|
|
1416 bam_fname = argv[optind++];
|
|
|
1417
|
|
|
1418 if ( !bam_fname )
|
|
|
1419 {
|
|
|
1420 if ( isatty(fileno((FILE *)stdin)) )
|
|
|
1421 error(NULL);
|
|
|
1422 bam_fname = "-";
|
|
|
1423 }
|
|
|
1424
|
|
|
1425 // Init structures
|
|
|
1426 // .. coverage bins and round buffer
|
|
|
1427 if ( stats->cov_step > stats->cov_max - stats->cov_min + 1 )
|
|
|
1428 {
|
|
|
1429 stats->cov_step = stats->cov_max - stats->cov_min;
|
|
|
1430 if ( stats->cov_step <= 0 )
|
|
|
1431 stats->cov_step = 1;
|
|
|
1432 }
|
|
|
1433 stats->ncov = 3 + (stats->cov_max-stats->cov_min) / stats->cov_step;
|
|
|
1434 stats->cov_max = stats->cov_min + ((stats->cov_max-stats->cov_min)/stats->cov_step +1)*stats->cov_step - 1;
|
|
|
1435 stats->cov = calloc(sizeof(uint64_t),stats->ncov);
|
|
|
1436 stats->cov_rbuf.size = stats->nbases*5;
|
|
|
1437 stats->cov_rbuf.buffer = calloc(sizeof(int32_t),stats->cov_rbuf.size);
|
|
|
1438 // .. bam
|
|
|
1439 if ((sam = samopen(bam_fname, in_mode, NULL)) == 0)
|
|
|
1440 error("Failed to open: %s\n", bam_fname);
|
|
|
1441 stats->sam = sam;
|
|
|
1442 if ( group_id ) init_group_id(stats, group_id);
|
|
|
1443 bam1_t *bam_line = bam_init1();
|
|
|
1444 // .. arrays
|
|
|
1445 stats->quals_1st = calloc(stats->nquals*stats->nbases,sizeof(uint64_t));
|
|
|
1446 stats->quals_2nd = calloc(stats->nquals*stats->nbases,sizeof(uint64_t));
|
|
|
1447 stats->gc_1st = calloc(stats->ngc,sizeof(uint64_t));
|
|
|
1448 stats->gc_2nd = calloc(stats->ngc,sizeof(uint64_t));
|
|
|
1449 stats->isize_inward = calloc(stats->nisize,sizeof(uint64_t));
|
|
|
1450 stats->isize_outward = calloc(stats->nisize,sizeof(uint64_t));
|
|
|
1451 stats->isize_other = calloc(stats->nisize,sizeof(uint64_t));
|
|
|
1452 stats->gcd = calloc(stats->ngcd,sizeof(gc_depth_t));
|
|
|
1453 stats->mpc_buf = stats->fai ? calloc(stats->nquals*stats->nbases,sizeof(uint64_t)) : NULL;
|
|
|
1454 stats->acgt_cycles = calloc(4*stats->nbases,sizeof(uint64_t));
|
|
|
1455 stats->read_lengths = calloc(stats->nbases,sizeof(uint64_t));
|
|
|
1456 stats->insertions = calloc(stats->nbases,sizeof(uint64_t));
|
|
|
1457 stats->deletions = calloc(stats->nbases,sizeof(uint64_t));
|
|
|
1458 stats->ins_cycles_1st = calloc(stats->nbases+1,sizeof(uint64_t));
|
|
|
1459 stats->ins_cycles_2nd = calloc(stats->nbases+1,sizeof(uint64_t));
|
|
|
1460 stats->del_cycles_1st = calloc(stats->nbases+1,sizeof(uint64_t));
|
|
|
1461 stats->del_cycles_2nd = calloc(stats->nbases+1,sizeof(uint64_t));
|
|
|
1462 realloc_rseq_buffer(stats);
|
|
|
1463 if ( targets )
|
|
|
1464 init_regions(stats, targets);
|
|
|
1465
|
|
|
1466 // Collect statistics
|
|
|
1467 if ( optind<argc )
|
|
|
1468 {
|
|
|
1469 // Collect stats in selected regions only
|
|
|
1470 bam_index_t *bam_idx = bam_index_load(bam_fname);
|
|
|
1471 if (bam_idx == 0)
|
|
|
1472 error("Random alignment retrieval only works for indexed BAM files.\n");
|
|
|
1473
|
|
|
1474 int i;
|
|
|
1475 for (i=optind; i<argc; i++)
|
|
|
1476 {
|
|
|
1477 int tid, beg, end;
|
|
|
1478 bam_parse_region(stats->sam->header, argv[i], &tid, &beg, &end);
|
|
|
1479 if ( tid < 0 ) continue;
|
|
|
1480 reset_regions(stats);
|
|
|
1481 bam_fetch(stats->sam->x.bam, bam_idx, tid, beg, end, stats, fetch_read);
|
|
|
1482 }
|
|
|
1483 bam_index_destroy(bam_idx);
|
|
|
1484 }
|
|
|
1485 else
|
|
|
1486 {
|
|
|
1487 // Stream through the entire BAM ignoring off-target regions if -t is given
|
|
|
1488 while (samread(sam,bam_line) >= 0)
|
|
|
1489 collect_stats(bam_line,stats);
|
|
|
1490 }
|
|
|
1491 round_buffer_flush(stats,-1);
|
|
|
1492
|
|
|
1493 output_stats(stats);
|
|
|
1494
|
|
|
1495 bam_destroy1(bam_line);
|
|
|
1496 samclose(stats->sam);
|
|
|
1497 if (stats->fai) fai_destroy(stats->fai);
|
|
|
1498 free(stats->cov_rbuf.buffer); free(stats->cov);
|
|
|
1499 free(stats->quals_1st); free(stats->quals_2nd);
|
|
|
1500 free(stats->gc_1st); free(stats->gc_2nd);
|
|
|
1501 free(stats->isize_inward); free(stats->isize_outward); free(stats->isize_other);
|
|
|
1502 free(stats->gcd);
|
|
|
1503 free(stats->rseq_buf);
|
|
|
1504 free(stats->mpc_buf);
|
|
|
1505 free(stats->acgt_cycles);
|
|
|
1506 free(stats->read_lengths);
|
|
|
1507 free(stats->insertions);
|
|
|
1508 free(stats->deletions);
|
|
|
1509 free(stats->ins_cycles_1st);
|
|
|
1510 free(stats->ins_cycles_2nd);
|
|
|
1511 free(stats->del_cycles_1st);
|
|
|
1512 free(stats->del_cycles_2nd);
|
|
|
1513 destroy_regions(stats);
|
|
|
1514 free(stats);
|
|
|
1515 if ( stats->rg_hash ) kh_destroy(kh_rg, stats->rg_hash);
|
|
|
1516
|
|
|
1517 return 0;
|
|
|
1518 }
|
|
|
1519
|
|
|
1520
|
|
|
1521
|
