|
0
|
1 #!/usr/bin/env perl
|
|
|
2 #
|
|
|
3 # Author: petr.danecek@sanger
|
|
|
4 #
|
|
|
5
|
|
|
6 use strict;
|
|
|
7 use warnings;
|
|
|
8 use Carp;
|
|
|
9
|
|
|
10 my $opts = parse_params();
|
|
|
11 parse_bamcheck($opts);
|
|
|
12 plot_qualities($opts);
|
|
|
13 plot_acgt_cycles($opts);
|
|
|
14 plot_gc($opts);
|
|
|
15 plot_gc_depth($opts);
|
|
|
16 plot_isize($opts);
|
|
|
17 plot_coverage($opts);
|
|
|
18 plot_mismatches_per_cycle($opts);
|
|
|
19 plot_indel_dist($opts);
|
|
|
20 plot_indel_cycles($opts);
|
|
|
21
|
|
|
22 exit;
|
|
|
23
|
|
|
24 #--------------------------------
|
|
|
25
|
|
|
26 sub error
|
|
|
27 {
|
|
|
28 my (@msg) = @_;
|
|
|
29 if ( scalar @msg ) { confess @msg; }
|
|
|
30 die
|
|
|
31 "Usage: plot-bamcheck [OPTIONS] file.bam.bc\n",
|
|
|
32 " plot-bamcheck -p outdir/ file.bam.bc\n",
|
|
|
33 "Options:\n",
|
|
|
34 " -k, --keep-files Do not remove temporary files.\n",
|
|
|
35 " -p, --prefix <path> The output files prefix, add a slash to create new directory.\n",
|
|
|
36 " -r, --ref-stats <file.fa.gc> Optional reference stats file with expected GC content (created with -s).\n",
|
|
|
37 " -s, --do-ref-stats <file.fa> Calculate reference sequence GC for later use with -r\n",
|
|
|
38 " -t, --targets <file.tab> Restrict -s to the listed regions (tab-delimited chr,from,to. 1-based, inclusive)\n",
|
|
|
39 " -h, -?, --help This help message.\n",
|
|
|
40 "\n";
|
|
|
41 }
|
|
|
42
|
|
|
43
|
|
|
44 sub parse_params
|
|
|
45 {
|
|
|
46 $0 =~ s{^.+/}{};
|
|
|
47 my $opts = { args=>join(' ',$0,@ARGV) };
|
|
|
48 while (defined(my $arg=shift(@ARGV)))
|
|
|
49 {
|
|
|
50 if ( $arg eq '-k' || $arg eq '--keep-files' ) { $$opts{keep_files}=1; next; }
|
|
|
51 if ( $arg eq '-r' || $arg eq '--ref-stats' ) { $$opts{ref_stats}=shift(@ARGV); next; }
|
|
|
52 if ( $arg eq '-s' || $arg eq '--do-ref-stats' ) { $$opts{do_ref_stats}=shift(@ARGV); next; }
|
|
|
53 if ( $arg eq '-t' || $arg eq '--targets' ) { $$opts{targets}=shift(@ARGV); next; }
|
|
|
54 if ( $arg eq '-p' || $arg eq '--prefix' ) { $$opts{prefix}=shift(@ARGV); next; }
|
|
|
55 if ( $arg eq '-?' || $arg eq '-h' || $arg eq '--help' ) { error(); }
|
|
|
56 if ( -e $arg ) { $$opts{bamcheck}=$arg; next; }
|
|
|
57 error("Unknown parameter or non-existent file \"$arg\". Run -h for help.\n");
|
|
|
58 }
|
|
|
59 if ( exists($$opts{do_ref_stats }) ) { do_ref_stats($opts); exit; }
|
|
|
60 if ( !exists($$opts{bamcheck}) ) { error("No bamcheck file?\n") }
|
|
|
61 if ( !exists($$opts{prefix}) ) { error("Expected -p parameter.\n") }
|
|
|
62 if ( $$opts{prefix}=~m{/$} ) { `mkdir -p $$opts{prefix}`; }
|
|
|
63 elsif ( !($$opts{prefix}=~/-$/) ) { $$opts{prefix} .= '-'; }
|
|
|
64 return $opts;
|
|
|
65 }
|
|
|
66
|
|
|
67
|
|
|
68 # Creates GC stats for either the whole reference or only on target regions for exome QC
|
|
|
69 sub do_ref_stats
|
|
|
70 {
|
|
|
71 my ($opts) = @_;
|
|
|
72
|
|
|
73
|
|
|
74 my %targets = ();
|
|
|
75 if ( exists($$opts{targets}) )
|
|
|
76 {
|
|
|
77 my ($prev_chr,$prev_pos);
|
|
|
78 open(my $fh,'<',$$opts{targets}) or error("$$opts{targets}: $!");
|
|
|
79 while (my $line=<$fh>)
|
|
|
80 {
|
|
|
81 if ( $line=~/^#/ ) { next; }
|
|
|
82 my ($chr,$from,$to) = split(/\s+/,$line);
|
|
|
83 chomp($to);
|
|
|
84 push @{$targets{$chr}}, $from,$to;
|
|
|
85 if ( !defined $prev_chr or $chr ne $prev_chr ) { $prev_chr=$chr; $prev_pos=$from }
|
|
|
86 if ( $prev_pos > $from ) { error("The file must be sorted: $$opts{targets}\n"); }
|
|
|
87 $prev_pos = $from;
|
|
|
88 }
|
|
|
89 close($fh);
|
|
|
90 }
|
|
|
91
|
|
|
92 my $_len = 60; # for now do only standard fasta's with 60 bases per line
|
|
|
93 my %gc_counts = ();
|
|
|
94 my ($skip_chr,$pos,$ireg,$regions);
|
|
|
95 open(my $fh,'<',$$opts{do_ref_stats}) or error("$$opts{do_ref_stats}: $!");
|
|
|
96 while (my $line=<$fh>)
|
|
|
97 {
|
|
|
98 if ( $line=~/^>/ )
|
|
|
99 {
|
|
|
100 if ( !scalar %targets ) { next; }
|
|
|
101
|
|
|
102 if ( !($line=~/>(\S+)/) ) { error("FIXME: could not determine chromosome name: $line"); }
|
|
|
103 if ( !exists($targets{$1}) ) { $skip_chr=$1; next; }
|
|
|
104 undef $skip_chr;
|
|
|
105 $pos = 0;
|
|
|
106 $ireg = 0;
|
|
|
107 $regions = $targets{$1};
|
|
|
108 }
|
|
|
109 if ( defined $skip_chr ) { next; }
|
|
|
110
|
|
|
111 # Only $_len sized lines are considered and no chopping for target regions.
|
|
|
112 chomp($line);
|
|
|
113 my $len = length($line);
|
|
|
114 if ( $len ne $_len ) { next; }
|
|
|
115
|
|
|
116 if ( scalar %targets )
|
|
|
117 {
|
|
|
118 while ( $ireg<@$regions && $$regions[$ireg+1]<=$pos ) { $ireg += 2; }
|
|
|
119 $pos += $len;
|
|
|
120 if ( $ireg==@$regions ) { next; }
|
|
|
121 if ( $pos < $$regions[$ireg] ) { next; }
|
|
|
122 }
|
|
|
123
|
|
|
124 my $gc_count = 0;
|
|
|
125 for (my $i=0; $i<$len; $i++)
|
|
|
126 {
|
|
|
127 my $base = substr($line,$i,1);
|
|
|
128 if ( $base eq 'g' || $base eq 'G' || $base eq 'c' || $base eq 'C' ) { $gc_count++; }
|
|
|
129 }
|
|
|
130 $gc_counts{$gc_count}++;
|
|
|
131 }
|
|
|
132
|
|
|
133 print "# Generated by $$opts{args}\n";
|
|
|
134 print "# The columns are: GC content bin, normalized frequency\n";
|
|
|
135 my $max;
|
|
|
136 for my $count (values %gc_counts)
|
|
|
137 {
|
|
|
138 if ( !defined $max or $count>$max ) { $max=$count; }
|
|
|
139 }
|
|
|
140 for my $gc (sort {$a<=>$b} keys %gc_counts)
|
|
|
141 {
|
|
|
142 if ( $gc==0 ) { next; }
|
|
|
143 printf "%f\t%f\n", $gc*100./$_len, $gc_counts{$gc}/$max;
|
|
|
144 }
|
|
|
145 }
|
|
|
146
|
|
|
147 sub plot
|
|
|
148 {
|
|
|
149 my ($cmdfile) = @_;
|
|
|
150 my $cmd = "gnuplot $cmdfile";
|
|
|
151 system($cmd);
|
|
|
152 if ( $? ) { error("The command exited with non-zero status $?:\n\t$cmd\n\n"); }
|
|
|
153 }
|
|
|
154
|
|
|
155
|
|
|
156 sub parse_bamcheck
|
|
|
157 {
|
|
|
158 my ($opts) = @_;
|
|
|
159 open(my $fh,'<',$$opts{bamcheck}) or error("$$opts{bamcheck}: $!");
|
|
|
160 my $line = <$fh>;
|
|
|
161 if ( !($line=~/^# This file was produced by bamcheck (\S+)/) ) { error("Sanity check failed: was this file generated by bamcheck?"); }
|
|
|
162 $$opts{dat}{version} = $1;
|
|
|
163 while ($line=<$fh>)
|
|
|
164 {
|
|
|
165 if ( $line=~/^#/ ) { next; }
|
|
|
166 my @items = split(/\t/,$line);
|
|
|
167 chomp($items[-1]);
|
|
|
168 if ( $items[0] eq 'SN' )
|
|
|
169 {
|
|
|
170 $$opts{dat}{$items[1]} = splice(@items,2);
|
|
|
171 next;
|
|
|
172 }
|
|
|
173 push @{$$opts{dat}{$items[0]}}, [splice(@items,1)];
|
|
|
174 }
|
|
|
175 close($fh);
|
|
|
176
|
|
|
177 # Check sanity
|
|
|
178 if ( !exists($$opts{dat}{'sequences:'}) or !$$opts{dat}{'sequences:'} )
|
|
|
179 {
|
|
|
180 error("Sanity check failed: no sequences found by bamcheck??\n");
|
|
|
181 }
|
|
|
182 }
|
|
|
183
|
|
|
184 sub older_than
|
|
|
185 {
|
|
|
186 my ($opts,$version) = @_;
|
|
|
187 my ($year,$month,$day) = split(/-/,$version);
|
|
|
188 $version = $$opts{dat}{version};
|
|
|
189 if ( !($version=~/\((\d+)-(\d+)-(\d+)\)$/) ) { return 1; }
|
|
|
190 if ( $1<$year ) { return 1; }
|
|
|
191 elsif ( $1>$year ) { return 0; }
|
|
|
192 if ( $2<$month ) { return 1; }
|
|
|
193 elsif ( $2>$month ) { return 0; }
|
|
|
194 if ( $3<$day ) { return 1; }
|
|
|
195 return 0;
|
|
|
196 }
|
|
|
197
|
|
|
198 sub get_defaults
|
|
|
199 {
|
|
|
200 my ($opts,$img_fname,%args) = @_;
|
|
|
201
|
|
|
202 if ( !($img_fname=~/\.png$/i) ) { error("FIXME: currently only PNG supported. (Easy to extend.)\n"); }
|
|
|
203
|
|
|
204 # Determine the gnuplot script file name
|
|
|
205 my $gp_file = $img_fname;
|
|
|
206 $gp_file =~ s{\.[^.]+$}{.gp};
|
|
|
207 if ( !($gp_file=~/.gp$/) ) { $gp_file .= '.gp'; }
|
|
|
208
|
|
|
209 # Determine the default title:
|
|
|
210 # 5446_6/5446_6.bam.bc.gp -> 5446_6
|
|
|
211 # test.aaa.png -> test.aaa
|
|
|
212 if ( !($$opts{bamcheck}=~m{([^/]+?)(?:\.bam)?(?:\.bc)?$}i) ) { error("FIXME: Could not determine the title from [$img_fname]\n"); }
|
|
|
213 my $title = $1;
|
|
|
214
|
|
|
215 my $dir = $gp_file;
|
|
|
216 $dir =~ s{/[^/]+$}{};
|
|
|
217 if ( $dir && $dir ne $gp_file ) { `mkdir -p $dir`; }
|
|
|
218
|
|
|
219 my $wh = exists($args{wh}) ? $args{wh} : '600,400';
|
|
|
220
|
|
|
221 open(my $fh,'>',$gp_file) or error("$gp_file: $!");
|
|
|
222 return {
|
|
|
223 title => $title,
|
|
|
224 gp => $gp_file,
|
|
|
225 img => $img_fname,
|
|
|
226 fh => $fh,
|
|
|
227 terminal => qq[set terminal png size $wh truecolor],
|
|
|
228 grid => 'set grid xtics ytics y2tics back lc rgb "#cccccc"',
|
|
|
229 };
|
|
|
230 }
|
|
|
231
|
|
|
232 sub percentile
|
|
|
233 {
|
|
|
234 my ($p,@vals) = @_;
|
|
|
235 my $N = 0;
|
|
|
236 for my $val (@vals) { $N += $val; }
|
|
|
237 my $n = $p*($N+1)/100.;
|
|
|
238 my $k = int($n);
|
|
|
239 my $d = $n-$k;
|
|
|
240 if ( $k<=0 ) { return 0; }
|
|
|
241 if ( $k>=$N ) { return scalar @vals-1; }
|
|
|
242 my $cnt;
|
|
|
243 for (my $i=0; $i<@vals; $i++)
|
|
|
244 {
|
|
|
245 $cnt += $vals[$i];
|
|
|
246 if ( $cnt>=$k ) { return $i; }
|
|
|
247 }
|
|
|
248 error("FIXME: this should not happen [percentile]\n");
|
|
|
249 }
|
|
|
250
|
|
|
251 sub plot_qualities
|
|
|
252 {
|
|
|
253 my ($opts) = @_;
|
|
|
254
|
|
|
255 if ( !exists($$opts{dat}{FFQ}) or !@{$$opts{dat}{FFQ}} ) { return; }
|
|
|
256
|
|
|
257 my $yrange = @{$$opts{dat}{FFQ}[0]} > 50 ? @{$$opts{dat}{FFQ}[0]} : 50;
|
|
|
258 my $is_paired = $$opts{dat}{'is paired:'};
|
|
|
259
|
|
|
260 # Average quality per cycle, forward and reverse reads in one plot
|
|
|
261 my $args = get_defaults($opts,"$$opts{prefix}quals.png");
|
|
|
262 my $fh = $$args{fh};
|
|
|
263 print $fh qq[
|
|
|
264 $$args{terminal}
|
|
|
265 set output "$$args{img}"
|
|
|
266 $$args{grid}
|
|
|
267 set ylabel "Average Quality"
|
|
|
268 set xlabel "Cycle"
|
|
|
269 set yrange [0:$yrange]
|
|
|
270 set title "$$args{title}"
|
|
|
271 plot '-' using 1:2 with lines title 'Forward reads' ] . ($is_paired ? q[, '-' using 1:2 with lines title 'Reverse reads'] : '') . q[
|
|
|
272 ];
|
|
|
273 my (@fp75,@fp50,@fmean);
|
|
|
274 my (@lp75,@lp50,@lmean);
|
|
|
275 my ($fmax,$fmax_qual,$fmax_cycle);
|
|
|
276 my ($lmax,$lmax_qual,$lmax_cycle);
|
|
|
277 for my $cycle (@{$$opts{dat}{FFQ}})
|
|
|
278 {
|
|
|
279 my $sum=0; my $n=0;
|
|
|
280 for (my $iqual=1; $iqual<@$cycle; $iqual++)
|
|
|
281 {
|
|
|
282 $sum += $$cycle[$iqual]*$iqual;
|
|
|
283 $n += $$cycle[$iqual];
|
|
|
284 if ( !defined $fmax or $fmax<$$cycle[$iqual] ) { $fmax=$$cycle[$iqual]; $fmax_qual=$iqual; $fmax_cycle=$$cycle[0]; }
|
|
|
285 }
|
|
|
286 my $p25 = percentile(25,(@$cycle)[1..$#$cycle]);
|
|
|
287 my $p50 = percentile(50,(@$cycle)[1..$#$cycle]);
|
|
|
288 my $p75 = percentile(75,(@$cycle)[1..$#$cycle]);
|
|
|
289 if ( !$n ) { next; }
|
|
|
290 push @fp75, "$$cycle[0]\t$p25\t$p75\n";
|
|
|
291 push @fp50, "$$cycle[0]\t$p50\n";
|
|
|
292 push @fmean, sprintf "%d\t%.2f\n", $$cycle[0],$sum/$n;
|
|
|
293 printf $fh $fmean[-1];
|
|
|
294 }
|
|
|
295 print $fh "end\n";
|
|
|
296 if ( $is_paired )
|
|
|
297 {
|
|
|
298 for my $cycle (@{$$opts{dat}{LFQ}})
|
|
|
299 {
|
|
|
300 my $sum=0; my $n=0;
|
|
|
301 for (my $iqual=1; $iqual<@$cycle; $iqual++)
|
|
|
302 {
|
|
|
303 $sum += $$cycle[$iqual]*$iqual;
|
|
|
304 $n += $$cycle[$iqual];
|
|
|
305 if ( !defined $lmax or $lmax<$$cycle[$iqual] ) { $lmax=$$cycle[$iqual]; $lmax_qual=$iqual; $lmax_cycle=$$cycle[0]; }
|
|
|
306 }
|
|
|
307 my $p25 = percentile(25,(@$cycle)[1..$#$cycle]);
|
|
|
308 my $p50 = percentile(50,(@$cycle)[1..$#$cycle]);
|
|
|
309 my $p75 = percentile(75,(@$cycle)[1..$#$cycle]);
|
|
|
310 if ( !$n ) { next; }
|
|
|
311 push @lp75, "$$cycle[0]\t$p25\t$p75\n";
|
|
|
312 push @lp50, "$$cycle[0]\t$p50\n";
|
|
|
313 push @lmean, sprintf "%d\t%.2f\n", $$cycle[0],$sum/$n;
|
|
|
314 printf $fh $lmean[-1];
|
|
|
315 }
|
|
|
316 print $fh "end\n";
|
|
|
317 }
|
|
|
318 close($fh);
|
|
|
319 plot($$args{gp});
|
|
|
320
|
|
|
321
|
|
|
322
|
|
|
323 # Average, mean and quality percentiles per cycle, forward and reverse reads in separate plots
|
|
|
324 $args = get_defaults($opts,"$$opts{prefix}quals2.png",wh=>'700,500');
|
|
|
325 $fh = $$args{fh};
|
|
|
326 print $fh qq[
|
|
|
327 $$args{terminal}
|
|
|
328 set output "$$args{img}"
|
|
|
329 $$args{grid}
|
|
|
330 set multiplot
|
|
|
331 set rmargin 0
|
|
|
332 set lmargin 0
|
|
|
333 set tmargin 0
|
|
|
334 set bmargin 0
|
|
|
335 set origin 0.1,0.1
|
|
|
336 set size 0.4,0.8
|
|
|
337 set yrange [0:$yrange]
|
|
|
338 set ylabel "Quality"
|
|
|
339 set xlabel "Cycle (fwd reads)"
|
|
|
340 plot '-' using 1:2:3 with filledcurve lt 1 lc rgb "#cccccc" t '25-75th percentile' , '-' using 1:2 with lines lc rgb "#000000" t 'Median', '-' using 1:2 with lines lt 1 t 'Mean'
|
|
|
341 ];
|
|
|
342 print $fh join('',@fp75),"end\n";
|
|
|
343 print $fh join('',@fp50),"end\n";
|
|
|
344 print $fh join('',@fmean),"end\n";
|
|
|
345 if ( $is_paired )
|
|
|
346 {
|
|
|
347 print $fh qq[
|
|
|
348 set origin 0.55,0.1
|
|
|
349 set size 0.4,0.8
|
|
|
350 unset ytics
|
|
|
351 set y2tics mirror
|
|
|
352 set yrange [0:$yrange]
|
|
|
353 unset ylabel
|
|
|
354 set xlabel "Cycle (rev reads)"
|
|
|
355 set label "$$args{title}" at screen 0.5,0.95 center
|
|
|
356 plot '-' using 1:2:3 with filledcurve lt 1 lc rgb "#cccccc" t '25-75th percentile' , '-' using 1:2 with lines lc rgb "#000000" t 'Median', '-' using 1:2 with lines lt 2 t 'Mean'
|
|
|
357 ];
|
|
|
358 print $fh join('',@lp75),"end\n";
|
|
|
359 print $fh join('',@lp50),"end\n";
|
|
|
360 print $fh join('',@lmean),"end\n";
|
|
|
361 }
|
|
|
362 close($fh);
|
|
|
363 plot($$args{gp});
|
|
|
364
|
|
|
365
|
|
|
366
|
|
|
367 # Quality distribution per cycle, the distribution is for each cycle plotted as a separate curve
|
|
|
368 $args = get_defaults($opts,"$$opts{prefix}quals3.png",wh=>'600,600');
|
|
|
369 $fh = $$args{fh};
|
|
|
370 my $nquals = @{$$opts{dat}{FFQ}[0]}-1;
|
|
|
371 my $ncycles = @{$$opts{dat}{FFQ}};
|
|
|
372 print $fh qq[
|
|
|
373 $$args{terminal}
|
|
|
374 set output "$$args{img}"
|
|
|
375 $$args{grid}
|
|
|
376 set multiplot
|
|
|
377 set rmargin 0
|
|
|
378 set lmargin 0
|
|
|
379 set tmargin 0
|
|
|
380 set bmargin 0
|
|
|
381 set origin 0.15,0.52
|
|
|
382 set size 0.8,0.4
|
|
|
383 set title "$$args{title}"
|
|
|
384 set ylabel "Frequency (fwd reads)"
|
|
|
385 set label "Cycle $fmax_cycle" at $fmax_qual+1,$fmax
|
|
|
386 unset xlabel
|
|
|
387 set xrange [0:$nquals]
|
|
|
388 set format x ""
|
|
|
389 ];
|
|
|
390 my @plots;
|
|
|
391 for (my $i=0; $i<$ncycles; $i++) { push @plots, q['-' using 1:2 with lines t ''] }
|
|
|
392 print $fh "plot ", join(",", @plots), "\n";
|
|
|
393 for my $cycle (@{$$opts{dat}{FFQ}})
|
|
|
394 {
|
|
|
395 for (my $iqual=1; $iqual<$nquals; $iqual++) { print $fh "$iqual\t$$cycle[$iqual]\n"; }
|
|
|
396 print $fh "end\n";
|
|
|
397 }
|
|
|
398 if ( $is_paired )
|
|
|
399 {
|
|
|
400 print $fh qq[
|
|
|
401 set origin 0.15,0.1
|
|
|
402 set size 0.8,0.4
|
|
|
403 unset title
|
|
|
404 unset format
|
|
|
405 set xtics
|
|
|
406 set xlabel "Quality"
|
|
|
407 unset label
|
|
|
408 set label "Cycle $lmax_cycle" at $lmax_qual+1,$lmax
|
|
|
409 set ylabel "Frequency (rev reads)"
|
|
|
410 ];
|
|
|
411 print $fh "plot ", join(",", @plots), "\n";
|
|
|
412 for my $cycle (@{$$opts{dat}{LFQ}})
|
|
|
413 {
|
|
|
414 for (my $iqual=1; $iqual<$nquals; $iqual++)
|
|
|
415 {
|
|
|
416 print $fh "$iqual\t$$cycle[$iqual]\n";
|
|
|
417 }
|
|
|
418 print $fh "end\n";
|
|
|
419 }
|
|
|
420 }
|
|
|
421 close($fh);
|
|
|
422 plot($$args{gp});
|
|
|
423
|
|
|
424
|
|
|
425 # Heatmap qualitites
|
|
|
426 $args = get_defaults($opts,"$$opts{prefix}quals-hm.png", wh=>'600,500');
|
|
|
427 $fh = $$args{fh};
|
|
|
428 my $max = defined $lmax && $lmax > $fmax ? $lmax : $fmax;
|
|
|
429 my @ytics;
|
|
|
430 for my $cycle (@{$$opts{dat}{FFQ}}) { if ( $$cycle[0]%10==0 ) { push @ytics,qq["$$cycle[0]" $$cycle[0]]; } }
|
|
|
431 my $ytics = join(',', @ytics);
|
|
|
432 print $fh qq[
|
|
|
433 $$args{terminal}
|
|
|
434 set output "$$args{img}"
|
|
|
435 unset key
|
|
|
436 unset colorbox
|
|
|
437 set palette defined (0 0 0 0, 1 0 0 1, 3 0 1 0, 4 1 0 0, 6 1 1 1)
|
|
|
438 set cbrange [0:$max]
|
|
|
439 set yrange [0:$ncycles]
|
|
|
440 set xrange [0:$nquals]
|
|
|
441 set view map
|
|
|
442 set multiplot
|
|
|
443 set rmargin 0
|
|
|
444 set lmargin 0
|
|
|
445 set tmargin 0
|
|
|
446 set bmargin 0
|
|
|
447 set origin 0,0.46
|
|
|
448 set size 0.95,0.6
|
|
|
449 set obj 1 rectangle behind from first 0,0 to first $nquals,$ncycles
|
|
|
450 set obj 1 fillstyle solid 1.0 fillcolor rgbcolor "black"
|
|
|
451 set ylabel "Cycle (fwd reads)" offset character -1,0
|
|
|
452 unset ytics
|
|
|
453 set ytics ($ytics)
|
|
|
454 unset xtics
|
|
|
455 set title "$$args{title}"
|
|
|
456 splot '-' matrix with image
|
|
|
457 ];
|
|
|
458 for my $cycle (@{$$opts{dat}{FFQ}})
|
|
|
459 {
|
|
|
460 for (my $iqual=1; $iqual<@$cycle; $iqual++) { print $fh "\t$$cycle[$iqual]"; }
|
|
|
461 print $fh "\n";
|
|
|
462 }
|
|
|
463 print $fh "end\nend\n";
|
|
|
464 @ytics = ();
|
|
|
465 for my $cycle (@{$$opts{dat}{LFQ}}) { if ( $$cycle[0]%10==0 ) { push @ytics,qq["$$cycle[0]" $$cycle[0]]; } }
|
|
|
466 $ytics = join(',', @ytics);
|
|
|
467 print $fh qq[
|
|
|
468 set origin 0,0.03
|
|
|
469 set size 0.95,0.6
|
|
|
470 set ylabel "Cycle (rev reads)" offset character -1,0
|
|
|
471 set xlabel "Base Quality"
|
|
|
472 unset title
|
|
|
473 unset ytics
|
|
|
474 set ytics ($ytics)
|
|
|
475 set xrange [0:$nquals]
|
|
|
476 set xtics
|
|
|
477 set colorbox vertical user origin first ($nquals+1),0 size screen 0.025,0.812
|
|
|
478 set cblabel "Number of bases"
|
|
|
479 splot '-' matrix with image
|
|
|
480 ];
|
|
|
481 for my $cycle (@{$$opts{dat}{LFQ}})
|
|
|
482 {
|
|
|
483 for (my $iqual=1; $iqual<@$cycle; $iqual++) { print $fh "\t$$cycle[$iqual]"; }
|
|
|
484 print $fh "\n";
|
|
|
485 }
|
|
|
486 print $fh "end\nend\n";
|
|
|
487 close($fh);
|
|
|
488 plot($$args{gp});
|
|
|
489 }
|
|
|
490
|
|
|
491
|
|
|
492 sub plot_acgt_cycles
|
|
|
493 {
|
|
|
494 my ($opts) = @_;
|
|
|
495
|
|
|
496 if ( !exists($$opts{dat}{GCC}) or !@{$$opts{dat}{GCC}} ) { return; }
|
|
|
497
|
|
|
498 my $args = get_defaults($opts,"$$opts{prefix}acgt-cycles.png");
|
|
|
499 my $fh = $$args{fh};
|
|
|
500 print $fh qq[
|
|
|
501 $$args{terminal}
|
|
|
502 set output "$$args{img}"
|
|
|
503 $$args{grid}
|
|
|
504 set style line 1 linecolor rgb "green"
|
|
|
505 set style line 2 linecolor rgb "red"
|
|
|
506 set style line 3 linecolor rgb "black"
|
|
|
507 set style line 4 linecolor rgb "blue"
|
|
|
508 set style increment user
|
|
|
509 set ylabel "Base content [%]"
|
|
|
510 set xlabel "Read Cycle"
|
|
|
511 set yrange [0:100]
|
|
|
512 set title "$$args{title}"
|
|
|
513 plot '-' w l ti 'A', '-' w l ti 'C', '-' w l ti 'G', '-' w l ti 'T'
|
|
|
514 ];
|
|
|
515 for my $base (1..4)
|
|
|
516 {
|
|
|
517 for my $cycle (@{$$opts{dat}{GCC}})
|
|
|
518 {
|
|
|
519 print $fh $$cycle[0]+1,"\t",$$cycle[$base],"\n";
|
|
|
520 }
|
|
|
521 print $fh "end\n";
|
|
|
522 }
|
|
|
523 close($fh);
|
|
|
524 plot($$args{gp});
|
|
|
525 }
|
|
|
526
|
|
|
527
|
|
|
528 sub plot_gc
|
|
|
529 {
|
|
|
530 my ($opts) = @_;
|
|
|
531
|
|
|
532 my $is_paired = $$opts{dat}{'is paired:'};
|
|
|
533 my $args = get_defaults($opts,"$$opts{prefix}gc-content.png");
|
|
|
534 my $fh = $$args{fh};
|
|
|
535 my ($gcl_max,$gcf_max,$lmax,$fmax);
|
|
|
536 for my $gc (@{$$opts{dat}{GCF}}) { if ( !defined $gcf_max or $gcf_max<$$gc[1] ) { $gcf_max=$$gc[1]; $fmax=$$gc[0]; } }
|
|
|
537 for my $gc (@{$$opts{dat}{GCL}}) { if ( !defined $gcl_max or $gcl_max<$$gc[1] ) { $gcl_max=$$gc[1]; $lmax=$$gc[0]; } }
|
|
|
538 my $gcmax = $is_paired && $gcl_max > $gcf_max ? $lmax : $fmax;
|
|
|
539 print $fh qq[
|
|
|
540 $$args{terminal}
|
|
|
541 set output "$$args{img}"
|
|
|
542 $$args{grid}
|
|
|
543 set title "$$args{title}"
|
|
|
544 set ylabel "Normalized Frequency"
|
|
|
545 set xlabel "GC Content [%]"
|
|
|
546 set yrange [0:1.1]
|
|
|
547 set label sprintf("%.1f",$gcmax) at $gcmax,1 front offset 1,0
|
|
|
548 plot ]
|
|
|
549 . (exists($$opts{ref_stats}) ? q['-' smooth csplines with lines lt 0 title 'Reference', ] : '')
|
|
|
550 . q['-' smooth csplines with lines lc 1 title 'First fragments' ]
|
|
|
551 . ($is_paired ? q[, '-' smooth csplines with lines lc 2 title 'Last fragments'] : '')
|
|
|
552 . q[
|
|
|
553 ];
|
|
|
554 if ( exists($$opts{ref_stats}) )
|
|
|
555 {
|
|
|
556 open(my $ref,'<',$$opts{ref_stats}) or error("$$opts{ref_stats}: $!");
|
|
|
557 while (my $line=<$ref>) { print $fh $line }
|
|
|
558 close($ref);
|
|
|
559 print $fh "end\n";
|
|
|
560 }
|
|
|
561 for my $cycle (@{$$opts{dat}{GCF}}) { printf $fh "%d\t%f\n", $$cycle[0],$$cycle[1]/$gcf_max; }
|
|
|
562 print $fh "end\n";
|
|
|
563 if ( $is_paired )
|
|
|
564 {
|
|
|
565 for my $cycle (@{$$opts{dat}{GCL}}) { printf $fh "%d\t%f\n", $$cycle[0],$$cycle[1]/$gcl_max; }
|
|
|
566 print $fh "end\n";
|
|
|
567 }
|
|
|
568 close($fh);
|
|
|
569 plot($$args{gp});
|
|
|
570 }
|
|
|
571
|
|
|
572
|
|
|
573 sub plot_gc_depth
|
|
|
574 {
|
|
|
575 my ($opts) = @_;
|
|
|
576
|
|
|
577 if ( !exists($$opts{dat}{GCD}) or !@{$$opts{dat}{GCD}} ) { return; }
|
|
|
578
|
|
|
579 # Find unique sequence percentiles for 30,40, and 50% GC content, just to draw x2tics.
|
|
|
580 my @tics = ( {gc=>30},{gc=>40},{gc=>50} );
|
|
|
581 for my $gc (@{$$opts{dat}{GCD}})
|
|
|
582 {
|
|
|
583 for my $tic (@tics)
|
|
|
584 {
|
|
|
585 my $diff = abs($$gc[0]-$$tic{gc});
|
|
|
586 if ( !exists($$tic{pr}) or $diff<$$tic{diff} ) { $$tic{pr}=$$gc[1]; $$tic{diff}=$diff; }
|
|
|
587 }
|
|
|
588 }
|
|
|
589
|
|
|
590 my @x2tics;
|
|
|
591 for my $tic (@tics) { push @x2tics, qq["$$tic{gc}" $$tic{pr}]; }
|
|
|
592 my $x2tics = join(',',@x2tics);
|
|
|
593
|
|
|
594 my $args = get_defaults($opts,"$$opts{prefix}gc-depth.png", wh=>'600,500');
|
|
|
595 my $fh = $$args{fh};
|
|
|
596 print $fh qq[
|
|
|
597 $$args{terminal}
|
|
|
598 set output "$$args{img}"
|
|
|
599 $$args{grid}
|
|
|
600 set ylabel "Mapped depth"
|
|
|
601 set xlabel "Percentile of mapped sequence ordered by GC content"
|
|
|
602 set x2label "GC Content [%]"
|
|
|
603 set title "$$args{title}"
|
|
|
604 set x2tics ($x2tics)
|
|
|
605 set xtics nomirror
|
|
|
606 set xrange [0.1:99.9]
|
|
|
607
|
|
|
608 plot '-' using 1:2:3 with filledcurve lt 1 lc rgb "#dedede" t '10-90th percentile' , \\
|
|
|
609 '-' using 1:2:3 with filledcurve lt 1 lc rgb "#bbdeff" t '25-75th percentile' , \\
|
|
|
610 '-' using 1:2 with lines lc rgb "#0084ff" t 'Median'
|
|
|
611 ];
|
|
|
612 for my $gc (@{$$opts{dat}{GCD}}) { print $fh "$$gc[1]\t$$gc[2]\t$$gc[6]\n"; } print $fh "end\n";
|
|
|
613 for my $gc (@{$$opts{dat}{GCD}}) { print $fh "$$gc[1]\t$$gc[3]\t$$gc[5]\n"; } print $fh "end\n";
|
|
|
614 for my $gc (@{$$opts{dat}{GCD}}) { print $fh "$$gc[1]\t$$gc[4]\n"; } print $fh "end\n";
|
|
|
615 close($fh);
|
|
|
616 plot($$args{gp});
|
|
|
617 }
|
|
|
618
|
|
|
619
|
|
|
620 sub plot_isize
|
|
|
621 {
|
|
|
622 my ($opts) = @_;
|
|
|
623
|
|
|
624 if ( !$$opts{dat}{'is paired:'} or !exists($$opts{dat}{IS}) or !@{$$opts{dat}{IS}} ) { return; }
|
|
|
625
|
|
|
626 my ($isize_max,$isize_cnt);
|
|
|
627 for my $isize (@{$$opts{dat}{IS}})
|
|
|
628 {
|
|
|
629 if ( !defined $isize_max or $isize_cnt<$$isize[1] ) { $isize_cnt=$$isize[1]; $isize_max=$$isize[0]; }
|
|
|
630 }
|
|
|
631
|
|
|
632 my $args = get_defaults($opts,"$$opts{prefix}insert-size.png");
|
|
|
633 my $fh = $$args{fh};
|
|
|
634 print $fh qq[
|
|
|
635 $$args{terminal}
|
|
|
636 set output "$$args{img}"
|
|
|
637 $$args{grid}
|
|
|
638 set rmargin 5
|
|
|
639 set label sprintf("%d",$isize_max) at $isize_max+10,$isize_cnt
|
|
|
640 set ylabel "Number of pairs"
|
|
|
641 set xlabel "Insert Size"
|
|
|
642 set title "$$args{title}"
|
|
|
643 plot \\
|
|
|
644 '-' with lines lc rgb 'black' title 'All pairs', \\
|
|
|
645 '-' with lines title 'Inward', \\
|
|
|
646 '-' with lines title 'Outward', \\
|
|
|
647 '-' with lines title 'Other'
|
|
|
648 ];
|
|
|
649 for my $isize (@{$$opts{dat}{IS}}) { print $fh "$$isize[0]\t$$isize[1]\n"; } print $fh "end\n";
|
|
|
650 for my $isize (@{$$opts{dat}{IS}}) { print $fh "$$isize[0]\t$$isize[2]\n"; } print $fh "end\n";
|
|
|
651 for my $isize (@{$$opts{dat}{IS}}) { print $fh "$$isize[0]\t$$isize[3]\n"; } print $fh "end\n";
|
|
|
652 for my $isize (@{$$opts{dat}{IS}}) { print $fh "$$isize[0]\t$$isize[4]\n"; } print $fh "end\n";
|
|
|
653 close($fh);
|
|
|
654 plot($$args{gp});
|
|
|
655 }
|
|
|
656
|
|
|
657
|
|
|
658 sub plot_coverage
|
|
|
659 {
|
|
|
660 my ($opts) = @_;
|
|
|
661
|
|
|
662 if ( !exists($$opts{dat}{COV}) or !@{$$opts{dat}{COV}} ) { return; }
|
|
|
663
|
|
|
664 my @vals;
|
|
|
665 for my $cov (@{$$opts{dat}{COV}}) { push @vals,$$cov[2]; }
|
|
|
666 my $i = percentile(99.8,@vals);
|
|
|
667 my $p99 = $$opts{dat}{COV}[$i][1];
|
|
|
668
|
|
|
669 my $args = get_defaults($opts,"$$opts{prefix}coverage.png");
|
|
|
670 my $fh = $$args{fh};
|
|
|
671 print $fh qq[
|
|
|
672 $$args{terminal}
|
|
|
673 set output "$$args{img}"
|
|
|
674 $$args{grid}
|
|
|
675 set ylabel "Number of mapped bases"
|
|
|
676 set xlabel "Coverage"
|
|
|
677 set style fill solid border -1
|
|
|
678 set title "$$args{title}"
|
|
|
679 set xrange [:$p99]
|
|
|
680 plot '-' with lines notitle
|
|
|
681 ];
|
|
|
682 for my $cov (@{$$opts{dat}{COV}})
|
|
|
683 {
|
|
|
684 if ( $$cov[2]==0 ) { next; }
|
|
|
685 print $fh "$$cov[1]\t$$cov[2]\n";
|
|
|
686 }
|
|
|
687 print $fh "end\n";
|
|
|
688 close($fh);
|
|
|
689 plot($$args{gp});
|
|
|
690 }
|
|
|
691
|
|
|
692
|
|
|
693 sub plot_mismatches_per_cycle
|
|
|
694 {
|
|
|
695 my ($opts) = @_;
|
|
|
696
|
|
|
697 if ( !exists($$opts{dat}{MPC}) or !@{$$opts{dat}{MPC}} ) { return; }
|
|
|
698 if ( older_than($opts,'2012-02-06') ) { plot_mismatches_per_cycle_old($opts); }
|
|
|
699
|
|
|
700 my $nquals = @{$$opts{dat}{MPC}[0]} - 2;
|
|
|
701 my $ncycles = @{$$opts{dat}{MPC}};
|
|
|
702 my ($style,$with);
|
|
|
703 if ( $ncycles>100 ) { $style = ''; $with = 'w l'; }
|
|
|
704 else { $style = 'set style data histogram; set style histogram rowstacked'; $with = ''; }
|
|
|
705
|
|
|
706 my $args = get_defaults($opts,"$$opts{prefix}mism-per-cycle.png");
|
|
|
707 my $fh = $$args{fh};
|
|
|
708 print $fh qq[
|
|
|
709 $$args{terminal}
|
|
|
710 set output "$$args{img}"
|
|
|
711 $$args{grid}
|
|
|
712 set style line 1 linecolor rgb "#e40000"
|
|
|
713 set style line 2 linecolor rgb "#ff9f00"
|
|
|
714 set style line 3 linecolor rgb "#eeee00"
|
|
|
715 set style line 4 linecolor rgb "#4ebd68"
|
|
|
716 set style line 5 linecolor rgb "#0061ff"
|
|
|
717 set style increment user
|
|
|
718 set key left top
|
|
|
719 $style
|
|
|
720 set ylabel "Number of mismatches"
|
|
|
721 set xlabel "Read Cycle"
|
|
|
722 set style fill solid border -1
|
|
|
723 set title "$$args{title}"
|
|
|
724 set xrange [-1:$ncycles]
|
|
|
725 plot '-' $with ti 'Base Quality>30', \\
|
|
|
726 '-' $with ti '30>=Q>20', \\
|
|
|
727 '-' $with ti '20>=Q>10', \\
|
|
|
728 '-' $with ti '10>=Q', \\
|
|
|
729 '-' $with ti "N's"
|
|
|
730 ];
|
|
|
731 for my $cycle (@{$$opts{dat}{MPC}})
|
|
|
732 {
|
|
|
733 my $sum; for my $idx (31..$#$cycle) { $sum += $$cycle[$idx]; }
|
|
|
734 print $fh "$sum\n";
|
|
|
735 }
|
|
|
736 print $fh "end\n";
|
|
|
737 for my $cycle (@{$$opts{dat}{MPC}})
|
|
|
738 {
|
|
|
739 my $sum; for my $idx (22..31) { $sum += $$cycle[$idx]; }
|
|
|
740 print $fh "$sum\n";
|
|
|
741 }
|
|
|
742 print $fh "end\n";
|
|
|
743 for my $cycle (@{$$opts{dat}{MPC}})
|
|
|
744 {
|
|
|
745 my $sum; for my $idx (12..21) { $sum += $$cycle[$idx]; }
|
|
|
746 print $fh "$sum\n";
|
|
|
747 }
|
|
|
748 print $fh "end\n";
|
|
|
749 for my $cycle (@{$$opts{dat}{MPC}})
|
|
|
750 {
|
|
|
751 my $sum; for my $idx (2..11) { $sum += $$cycle[$idx]; }
|
|
|
752 print $fh "$sum\n";
|
|
|
753 }
|
|
|
754 print $fh "end\n";
|
|
|
755 for my $cycle (@{$$opts{dat}{MPC}}) { print $fh "$$cycle[1]\n"; }
|
|
|
756 print $fh "end\n";
|
|
|
757 close($fh);
|
|
|
758 plot($$args{gp});
|
|
|
759 }
|
|
|
760
|
|
|
761 sub plot_indel_dist
|
|
|
762 {
|
|
|
763 my ($opts) = @_;
|
|
|
764
|
|
|
765 if ( !exists($$opts{dat}{ID}) or !@{$$opts{dat}{ID}} ) { return; }
|
|
|
766
|
|
|
767 my $args = get_defaults($opts,"$$opts{prefix}indel-dist.png");
|
|
|
768 my $fh = $$args{fh};
|
|
|
769 print $fh qq[
|
|
|
770 $$args{terminal}
|
|
|
771 set output "$$args{img}"
|
|
|
772 $$args{grid}
|
|
|
773 set style line 1 linetype 1 linecolor rgb "red"
|
|
|
774 set style line 2 linetype 2 linecolor rgb "black"
|
|
|
775 set style line 3 linetype 3 linecolor rgb "green"
|
|
|
776 set style increment user
|
|
|
777 set ylabel "Indel count [log]"
|
|
|
778 set xlabel "Indel length"
|
|
|
779 set y2label "Insertions/Deletions ratio"
|
|
|
780 set log y
|
|
|
781 set y2tics nomirror
|
|
|
782 set ytics nomirror
|
|
|
783 set title "$$args{title}"
|
|
|
784 plot '-' w l ti 'Insertions', '-' w l ti 'Deletions', '-' axes x1y2 w l ti "Ins/Dels ratio"
|
|
|
785 ];
|
|
|
786 for my $len (@{$$opts{dat}{ID}}) { print $fh "$$len[0]\t$$len[1]\n"; } print $fh "end\n";
|
|
|
787 for my $len (@{$$opts{dat}{ID}}) { print $fh "$$len[0]\t$$len[2]\n"; } print $fh "end\n";
|
|
|
788 for my $len (@{$$opts{dat}{ID}}) { printf $fh "%d\t%f\n", $$len[0],$$len[2]?$$len[1]/$$len[2]:0; } print $fh "end\n";
|
|
|
789 close($fh);
|
|
|
790 plot($$args{gp});
|
|
|
791 }
|
|
|
792
|
|
|
793 sub plot_indel_cycles
|
|
|
794 {
|
|
|
795 my ($opts) = @_;
|
|
|
796
|
|
|
797 if ( !exists($$opts{dat}{IC}) or !@{$$opts{dat}{IC}} ) { return; }
|
|
|
798
|
|
|
799 my $args = get_defaults($opts,"$$opts{prefix}indel-cycles.png");
|
|
|
800 my $fh = $$args{fh};
|
|
|
801 print $fh qq[
|
|
|
802 $$args{terminal}
|
|
|
803 set output "$$args{img}"
|
|
|
804 $$args{grid}
|
|
|
805 set style line 1 linetype 1 linecolor rgb "red"
|
|
|
806 set style line 2 linetype 2 linecolor rgb "black"
|
|
|
807 set style line 3 linetype 3 linecolor rgb "green"
|
|
|
808 set style line 4 linetype 4 linecolor rgb "blue"
|
|
|
809 set style increment user
|
|
|
810 set ylabel "Indel count"
|
|
|
811 set xlabel "Read Cycle"
|
|
|
812 set title "$$args{title}"
|
|
|
813 plot '-' w l ti 'Insertions (fwd)', '' w l ti 'Insertions (rev)', '' w l ti 'Deletions (fwd)', '' w l ti 'Deletions (rev)'
|
|
|
814 ];
|
|
|
815 for my $len (@{$$opts{dat}{IC}}) { print $fh "$$len[0]\t$$len[1]\n"; } print $fh "end\n";
|
|
|
816 for my $len (@{$$opts{dat}{IC}}) { print $fh "$$len[0]\t$$len[2]\n"; } print $fh "end\n";
|
|
|
817 for my $len (@{$$opts{dat}{IC}}) { print $fh "$$len[0]\t$$len[3]\n"; } print $fh "end\n";
|
|
|
818 for my $len (@{$$opts{dat}{IC}}) { print $fh "$$len[0]\t$$len[4]\n"; } print $fh "end\n";
|
|
|
819 close($fh);
|
|
|
820 plot($$args{gp});
|
|
|
821 }
|
|
|
822
|
|
|
823
|
|
|
824
|
|
|
825
|
|
|
826
|
|
|
827
|
|
|
828
|
|
|
829 sub has_values
|
|
|
830 {
|
|
|
831 my ($opts,@tags) = @_;
|
|
|
832 for my $tag (@tags)
|
|
|
833 {
|
|
|
834 my (@lines) = `cat $$opts{bamcheck} | grep ^$tag | wc -l`;
|
|
|
835 chomp($lines[0]);
|
|
|
836 if ( $lines[0]<2 ) { return 0; }
|
|
|
837 }
|
|
|
838 return 1;
|
|
|
839 }
|
|
|
840
|
|
|
841 sub plot_mismatches_per_cycle_old
|
|
|
842 {
|
|
|
843 my ($opts) = @_;
|
|
|
844
|
|
|
845 my $args = get_defaults($opts,"$$opts{prefix}mism-per-cycle.png");
|
|
|
846 my ($nquals) = `grep ^MPC $$opts{bamcheck} | awk '\$2==1' | sed 's,\\t,\\n,g' | wc -l`;
|
|
|
847 my ($ncycles) = `grep ^MPC $$opts{bamcheck} | wc -l`;
|
|
|
848 chomp($nquals);
|
|
|
849 chomp($ncycles);
|
|
|
850 $nquals--;
|
|
|
851 $ncycles--;
|
|
|
852 my @gr0_15 = (2..17);
|
|
|
853 my @gr16_30 = (18..32);
|
|
|
854 my @gr31_n = (33..$nquals);
|
|
|
855 my $gr0_15 = '$'. join('+$',@gr0_15);
|
|
|
856 my $gr16_30 = '$'. join('+$',@gr16_30);
|
|
|
857 my $gr31_n = '$'. join('+$',@gr31_n);
|
|
|
858
|
|
|
859 open(my $fh,'>',$$args{gp}) or error("$$args{gp}: $!");
|
|
|
860 print $fh q[
|
|
|
861 set terminal png size 600,400 truecolor font "DejaVuSansMono,9"
|
|
|
862 set output "] . $$args{img} . q["
|
|
|
863
|
|
|
864 set key left top
|
|
|
865 set style data histogram
|
|
|
866 set style histogram rowstacked
|
|
|
867
|
|
|
868 set grid back lc rgb "#aaaaaa"
|
|
|
869 set ylabel "Number of mismatches"
|
|
|
870 set xlabel "Read Cycle"
|
|
|
871 set style fill solid border -1
|
|
|
872 set title "] . $$args{title} . qq["
|
|
|
873 set xrange [-1:$ncycles]
|
|
|
874
|
|
|
875 plot '< grep ^MPC $$opts{bamcheck} | cut -f 2-' using ($gr31_n) ti 'Base Quality>30', '' using ($gr16_30) ti '30>=Q>15', '' using ($gr0_15) ti '15>=Q'
|
|
|
876 ];
|
|
|
877 close($fh);
|
|
|
878
|
|
|
879 plot($$args{gp});
|
|
|
880 }
|
|
|
881
|
|
|
882
|
