misc/plot-bamcheck

   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                 unset ylabel
 353                 set xlabel "Cycle (rev reads)"
 354                 set label "$$args{title}" at screen 0.5,0.95 center
 355                 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'
 356             ];
 357         print $fh join('',@lp75),"end\n";
 358         print $fh join('',@lp50),"end\n";
 359         print $fh join('',@lmean),"end\n";
 360     }
 361     close($fh);
 362     plot($$args{gp});
 363
 364
 365
 366     # Quality distribution per cycle, the distribution is for each cycle plotted as a separate curve
 367     $args = get_defaults($opts,"$$opts{prefix}quals3.png",wh=>'600,600');
 368     $fh = $$args{fh};
 369     my $nquals = @{$$opts{dat}{FFQ}[0]}-1;
 370     my $ncycles = @{$$opts{dat}{FFQ}};
 371     print $fh qq[
 372             $$args{terminal}
 373             set output "$$args{img}"
 374             $$args{grid}
 375             set multiplot
 376             set rmargin 0
 377             set lmargin 0
 378             set tmargin 0
 379             set bmargin 0
 380             set origin 0.15,0.52
 381             set size 0.8,0.4
 382             set title "$$args{title}"
 383             set ylabel "Frequency (fwd reads)"
 384             set label "Cycle $fmax_cycle" at $fmax_qual+1,$fmax
 385             unset xlabel
 386             set xrange [0:$nquals]
 387             set format x ""
 388         ];
 389     my @plots;
 390     for (my $i=0; $i<$ncycles; $i++) { push @plots, q['-' using 1:2 with lines t ''] }
 391     print $fh "plot ", join(",", @plots), "\n";
 392     for my $cycle (@{$$opts{dat}{FFQ}})
 393     {
 394         for (my $iqual=1; $iqual<$nquals; $iqual++) { print $fh "$iqual\t$$cycle[$iqual]\n"; }
 395         print $fh "end\n";
 396     }
 397     if ( $is_paired )
 398     {
 399         print $fh qq[
 400                 set origin 0.15,0.1
 401                 set size 0.8,0.4
 402                 unset title
 403                 unset format
 404                 set xtics
 405                 set xlabel "Quality"
 406                 unset label
 407                 set label "Cycle $lmax_cycle" at $lmax_qual+1,$lmax
 408                 set ylabel "Frequency (rev reads)"
 409             ];
 410         print $fh "plot ", join(",", @plots), "\n";
 411         for my $cycle (@{$$opts{dat}{LFQ}})
 412         {
 413             for (my $iqual=1; $iqual<$nquals; $iqual++)
 414             {
 415                 print $fh "$iqual\t$$cycle[$iqual]\n";
 416             }
 417             print $fh "end\n";
 418         }
 419     }
 420     close($fh);
 421     plot($$args{gp});
 422
 423
 424     # Heatmap qualitites
 425     $args = get_defaults($opts,"$$opts{prefix}quals-hm.png", wh=>'600,500');
 426     $fh = $$args{fh};
 427     my $max = defined $lmax && $lmax > $fmax ? $lmax : $fmax;
 428     my @ytics;
 429     for my $cycle (@{$$opts{dat}{FFQ}}) { if ( $$cycle[0]%10==0 ) { push @ytics,qq["$$cycle[0]" $$cycle[0]]; } }
 430     my $ytics = join(',', @ytics);
 431     print $fh qq[
 432             $$args{terminal}
 433             set output "$$args{img}"
 434             unset key
 435             unset colorbox
 436             set palette defined (0 0 0 0, 1 0 0 1, 3 0 1 0, 4 1 0 0, 6 1 1 1)
 437             set cbrange [0:$max]
 438             set yrange  [0:$ncycles]
 439             set xrange  [0:$nquals]
 440             set view map
 441             set multiplot
 442             set rmargin 0
 443             set lmargin 0
 444             set tmargin 0
 445             set bmargin 0
 446             set origin 0,0.46
 447             set size 0.95,0.6
 448             set obj 1 rectangle behind from first 0,0 to first $nquals,$ncycles
 449             set obj 1 fillstyle solid 1.0 fillcolor rgbcolor "black"
 450             set ylabel "Cycle (fwd reads)" offset character -1,0
 451             unset ytics
 452             set ytics ($ytics)
 453             unset xtics
 454             set title "$$args{title}"
 455             splot '-' matrix with image
 456         ];
 457     for my $cycle (@{$$opts{dat}{FFQ}})
 458     {
 459         for (my $iqual=1; $iqual<@$cycle; $iqual++) { print $fh "\t$$cycle[$iqual]"; }
 460         print $fh "\n";
 461     }
 462     print $fh "end\nend\n";
 463     @ytics = ();
 464     for my $cycle (@{$$opts{dat}{LFQ}}) { if ( $$cycle[0]%10==0 ) { push @ytics,qq["$$cycle[0]" $$cycle[0]]; } }
 465     $ytics = join(',', @ytics);
 466     print $fh qq[
 467             set origin 0,0.03
 468             set size 0.95,0.6
 469             set ylabel "Cycle (rev reads)" offset character -1,0
 470             set xlabel "Base Quality"
 471             unset title
 472             unset ytics
 473             set ytics ($ytics)
 474             set xrange  [0:$nquals]
 475             set xtics
 476             set colorbox vertical user origin first ($nquals+1),0 size screen 0.025,0.812
 477             set cblabel "Number of bases"
 478             splot '-' matrix with image
 479         ];
 480     for my $cycle (@{$$opts{dat}{LFQ}})
 481     {
 482         for (my $iqual=1; $iqual<@$cycle; $iqual++) { print $fh "\t$$cycle[$iqual]"; }
 483         print $fh "\n";
 484     }
 485     print $fh "end\nend\n";
 486     close($fh);
 487     plot($$args{gp});
 488 }
 489
 490
 491 sub plot_acgt_cycles
 492 {
 493     my ($opts) = @_;
 494
 495     if ( !exists($$opts{dat}{GCC}) or !@{$$opts{dat}{GCC}} ) { return; }
 496
 497     my $args = get_defaults($opts,"$$opts{prefix}acgt-cycles.png");
 498     my $fh = $$args{fh};
 499     print $fh qq[
 500             $$args{terminal}
 501             set output "$$args{img}"
 502             $$args{grid}
 503             set style line 1 linecolor rgb "green"
 504             set style line 2 linecolor rgb "red"
 505             set style line 3 linecolor rgb "black"
 506             set style line 4 linecolor rgb "blue"
 507             set style increment user
 508             set ylabel "Base content [%]"
 509             set xlabel "Read Cycle"
 510             set yrange [0:100]
 511             set title "$$args{title}"
 512             plot '-' w l ti 'A', '-' w l ti 'C', '-' w l ti 'G', '-' w l ti 'T'
 513         ];
 514     for my $base (1..4)
 515     {
 516         for my $cycle (@{$$opts{dat}{GCC}})
 517         {
 518             print $fh $$cycle[0]+1,"\t",$$cycle[$base],"\n";
 519         }
 520         print $fh "end\n";
 521     }
 522     close($fh);
 523     plot($$args{gp});
 524 }
 525
 526
 527 sub plot_gc
 528 {
 529     my ($opts) = @_;
 530
 531     my $is_paired = $$opts{dat}{'is paired:'};
 532     my $args = get_defaults($opts,"$$opts{prefix}gc-content.png");
 533     my $fh = $$args{fh};
 534     my ($gcl_max,$gcf_max,$lmax,$fmax);
 535     for my $gc (@{$$opts{dat}{GCF}}) { if ( !defined $gcf_max or $gcf_max<$$gc[1] ) { $gcf_max=$$gc[1]; $fmax=$$gc[0]; } }
 536     for my $gc (@{$$opts{dat}{GCL}}) { if ( !defined $gcl_max or $gcl_max<$$gc[1] ) { $gcl_max=$$gc[1]; $lmax=$$gc[0]; } }
 537     my $gcmax = $is_paired && $gcl_max > $gcf_max ? $lmax : $fmax;
 538     print $fh qq[
 539             $$args{terminal}
 540             set output "$$args{img}"
 541             $$args{grid}
 542             set title "$$args{title}"
 543             set ylabel "Normalized Frequency"
 544             set xlabel "GC Content [%]"
 545             set yrange [0:1.1]
 546             set label sprintf("%.1f",$gcmax) at $gcmax,1 front offset 1,0
 547             plot ]
 548                  . (exists($$opts{ref_stats}) ? q['-' smooth csplines with lines lt 0 title 'Reference', ] : '')
 549                  . q['-' smooth csplines with lines lc 1 title 'First fragments' ]
 550                  . ($is_paired ? q[, '-' smooth csplines with lines lc 2 title 'Last fragments'] : '')
 551                  . q[
 552         ];
 553     if ( exists($$opts{ref_stats}) )
 554     {
 555         open(my $ref,'<',$$opts{ref_stats}) or error("$$opts{ref_stats}: $!");
 556         while (my $line=<$ref>) { print $fh $line }
 557         close($ref);
 558         print $fh "end\n";
 559     }
 560     for my $cycle (@{$$opts{dat}{GCF}}) { printf $fh "%d\t%f\n", $$cycle[0],$$cycle[1]/$gcf_max; }
 561     print $fh "end\n";
 562     if ( $is_paired )
 563     {
 564         for my $cycle (@{$$opts{dat}{GCL}}) { printf $fh "%d\t%f\n", $$cycle[0],$$cycle[1]/$gcl_max; }
 565         print $fh "end\n";
 566     }
 567     close($fh);
 568     plot($$args{gp});
 569 }
 570
 571
 572 sub plot_gc_depth
 573 {
 574     my ($opts) = @_;
 575
 576     if ( !exists($$opts{dat}{GCD}) or !@{$$opts{dat}{GCD}} ) { return; }
 577
 578     # Find unique sequence percentiles for 30,40, and 50% GC content, just to draw x2tics.
 579     my @tics = ( {gc=>30},{gc=>40},{gc=>50} );
 580     for my $gc (@{$$opts{dat}{GCD}})
 581     {
 582         for my $tic (@tics)
 583         {
 584             my $diff = abs($$gc[0]-$$tic{gc});
 585             if ( !exists($$tic{pr}) or $diff<$$tic{diff} ) { $$tic{pr}=$$gc[1]; $$tic{diff}=$diff; }
 586         }
 587     }
 588
 589     my @x2tics;
 590     for my $tic (@tics) { push @x2tics, qq["$$tic{gc}" $$tic{pr}]; }
 591     my $x2tics = join(',',@x2tics);
 592
 593     my $args = get_defaults($opts,"$$opts{prefix}gc-depth.png", wh=>'600,500');
 594     my $fh = $$args{fh};
 595     print $fh qq[
 596             $$args{terminal}
 597             set output "$$args{img}"
 598             $$args{grid}
 599             set ylabel "Mapped depth"
 600             set xlabel "Percentile of mapped sequence ordered by GC content"
 601             set x2label "GC Content [%]"
 602             set title "$$args{title}"
 603             set x2tics ($x2tics)
 604             set xtics nomirror
 605             set xrange [0.1:99.9]
 606
 607             plot '-' using 1:2:3 with filledcurve lt 1 lc rgb "#dedede" t '10-90th percentile' , \\
 608                  '-' using 1:2:3 with filledcurve lt 1 lc rgb "#bbdeff" t '25-75th percentile' , \\
 609                  '-' using 1:2 with lines lc rgb "#0084ff" t 'Median'
 610         ];
 611     for my $gc (@{$$opts{dat}{GCD}}) { print $fh "$$gc[1]\t$$gc[2]\t$$gc[6]\n"; } print $fh "end\n";
 612     for my $gc (@{$$opts{dat}{GCD}}) { print $fh "$$gc[1]\t$$gc[3]\t$$gc[5]\n"; } print $fh "end\n";
 613     for my $gc (@{$$opts{dat}{GCD}}) { print $fh "$$gc[1]\t$$gc[4]\n"; } print $fh "end\n";
 614     close($fh);
 615     plot($$args{gp});
 616 }
 617
 618
 619 sub plot_isize
 620 {
 621     my ($opts) = @_;
 622
 623     if ( !$$opts{dat}{'is paired:'} or !exists($$opts{dat}{IS}) or !@{$$opts{dat}{IS}} ) { return; }
 624
 625     my ($isize_max,$isize_cnt);
 626     for my $isize (@{$$opts{dat}{IS}})
 627     {
 628         if ( !defined $isize_max or $isize_cnt<$$isize[1] ) { $isize_cnt=$$isize[1]; $isize_max=$$isize[0]; }
 629     }
 630
 631     my $args = get_defaults($opts,"$$opts{prefix}insert-size.png");
 632     my $fh = $$args{fh};
 633     print $fh qq[
 634             $$args{terminal}
 635             set output "$$args{img}"
 636             $$args{grid}
 637             set rmargin 5
 638             set label sprintf("%d",$isize_max) at $isize_max+10,$isize_cnt
 639             set ylabel  "Number of pairs"
 640             set xlabel  "Insert Size"
 641             set title "$$args{title}"
 642             plot \\
 643                 '-' with lines lc rgb 'black' title 'All pairs', \\
 644                 '-' with lines title 'Inward', \\
 645                 '-' with lines title 'Outward', \\
 646                 '-' with lines title 'Other'
 647         ];
 648     for my $isize (@{$$opts{dat}{IS}}) { print $fh "$$isize[0]\t$$isize[1]\n"; } print $fh "end\n";
 649     for my $isize (@{$$opts{dat}{IS}}) { print $fh "$$isize[0]\t$$isize[2]\n"; } print $fh "end\n";
 650     for my $isize (@{$$opts{dat}{IS}}) { print $fh "$$isize[0]\t$$isize[3]\n"; } print $fh "end\n";
 651     for my $isize (@{$$opts{dat}{IS}}) { print $fh "$$isize[0]\t$$isize[4]\n"; } print $fh "end\n";
 652     close($fh);
 653     plot($$args{gp});
 654 }
 655
 656
 657 sub plot_coverage
 658 {
 659     my ($opts) = @_;
 660
 661     if ( !exists($$opts{dat}{COV}) or !@{$$opts{dat}{COV}} ) { return; }
 662
 663     my @vals;
 664     for my $cov (@{$$opts{dat}{COV}}) { push @vals,$$cov[2]; }
 665     my $p99 = percentile(99.8,@vals);
 666
 667     my $args = get_defaults($opts,"$$opts{prefix}coverage.png");
 668     my $fh = $$args{fh};
 669     print $fh qq[
 670             $$args{terminal}
 671             set output "$$args{img}"
 672             $$args{grid}
 673             set ylabel "Number of mapped bases"
 674             set xlabel "Coverage"
 675             set style fill solid border -1
 676             set title "$$args{title}"
 677             set xrange [:$p99]
 678             plot '-' with lines notitle
 679         ];
 680     for my $cov (@{$$opts{dat}{COV}})
 681     {
 682         if ( $$cov[2]==0 ) { next; }
 683         print $fh "$$cov[1]\t$$cov[2]\n";
 684     }
 685     print $fh "end\n";
 686     close($fh);
 687     plot($$args{gp});
 688 }
 689
 690
 691 sub plot_mismatches_per_cycle
 692 {
 693     my ($opts) = @_;
 694
 695     if ( !exists($$opts{dat}{MPC}) or !@{$$opts{dat}{MPC}} ) { return; }
 696     if ( older_than($opts,'2012-02-06') ) { plot_mismatches_per_cycle_old($opts); }
 697
 698     my $nquals = @{$$opts{dat}{MPC}[0]} - 2;
 699     my $ncycles = @{$$opts{dat}{MPC}};
 700     my ($style,$with);
 701     if ( $ncycles>100 ) { $style = ''; $with = 'w l'; }
 702     else { $style = 'set style data histogram; set style histogram rowstacked'; $with = ''; }
 703
 704     my $args = get_defaults($opts,"$$opts{prefix}mism-per-cycle.png");
 705     my $fh = $$args{fh};
 706     print $fh qq[
 707             $$args{terminal}
 708             set output "$$args{img}"
 709             $$args{grid}
 710             set style line 1 linecolor rgb "#e40000"
 711             set style line 2 linecolor rgb "#ff9f00"
 712             set style line 3 linecolor rgb "#eeee00"
 713             set style line 4 linecolor rgb "#4ebd68"
 714             set style line 5 linecolor rgb "#0061ff"
 715             set style increment user
 716             set key left top
 717             $style
 718             set ylabel "Number of mismatches"
 719             set xlabel "Read Cycle"
 720             set style fill solid border -1
 721             set title "$$args{title}"
 722             set xrange [-1:$ncycles]
 723             plot '-' $with ti 'Base Quality>30', \\
 724                  '-' $with ti '30>=Q>20', \\
 725                  '-' $with ti '20>=Q>10', \\
 726                  '-' $with ti '10>=Q', \\
 727                  '-' $with ti "N's"
 728         ];
 729     for my $cycle (@{$$opts{dat}{MPC}})
 730     {
 731         my $sum; for my $idx (31..$#$cycle) { $sum += $$cycle[$idx]; }
 732         print $fh "$sum\n";
 733     }
 734     print $fh "end\n";
 735     for my $cycle (@{$$opts{dat}{MPC}})
 736     {
 737         my $sum; for my $idx (22..31) { $sum += $$cycle[$idx]; }
 738         print $fh "$sum\n";
 739     }
 740     print $fh "end\n";
 741     for my $cycle (@{$$opts{dat}{MPC}})
 742     {
 743         my $sum; for my $idx (12..21) { $sum += $$cycle[$idx]; }
 744         print $fh "$sum\n";
 745     }
 746     print $fh "end\n";
 747     for my $cycle (@{$$opts{dat}{MPC}})
 748     {
 749         my $sum; for my $idx (2..11) { $sum += $$cycle[$idx]; }
 750         print $fh "$sum\n";
 751     }
 752     print $fh "end\n";
 753     for my $cycle (@{$$opts{dat}{MPC}}) { print $fh "$$cycle[1]\n"; }
 754     print $fh "end\n";
 755     close($fh);
 756     plot($$args{gp});
 757 }
 758
 759 sub plot_indel_dist
 760 {
 761     my ($opts) = @_;
 762
 763     if ( !exists($$opts{dat}{ID}) or !@{$$opts{dat}{ID}} ) { return; }
 764
 765     my $args = get_defaults($opts,"$$opts{prefix}indel-dist.png");
 766     my $fh = $$args{fh};
 767     print $fh qq[
 768         $$args{terminal}
 769         set output "$$args{img}"
 770         $$args{grid}
 771         set style line 1 linetype 1  linecolor rgb "red"
 772         set style line 2 linetype 2  linecolor rgb "black"
 773         set style line 3 linetype 3  linecolor rgb "green"
 774         set style increment user
 775         set ylabel "Indel count [log]"
 776         set xlabel "Indel length"
 777         set y2label "Insertions/Deletions ratio"
 778         set log y
 779         set y2tics nomirror
 780         set ytics nomirror
 781         set title "$$args{title}"
 782         plot '-' w l ti 'Insertions', '-' w l ti 'Deletions', '-' axes x1y2 w l ti "Ins/Dels ratio"
 783     ];
 784     for my $len (@{$$opts{dat}{ID}}) { print $fh "$$len[0]\t$$len[1]\n"; } print $fh "end\n";
 785     for my $len (@{$$opts{dat}{ID}}) { print $fh "$$len[0]\t$$len[2]\n"; } print $fh "end\n";
 786     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";
 787     close($fh);
 788     plot($$args{gp});
 789 }
 790
 791 sub plot_indel_cycles
 792 {
 793     my ($opts) = @_;
 794
 795     if ( !exists($$opts{dat}{IC}) or !@{$$opts{dat}{IC}} ) { return; }
 796
 797     my $args = get_defaults($opts,"$$opts{prefix}indel-cycles.png");
 798     my $fh = $$args{fh};
 799     print $fh qq[
 800         $$args{terminal}
 801         set output "$$args{img}"
 802         $$args{grid}
 803         set style line 1 linetype 1  linecolor rgb "red"
 804         set style line 2 linetype 2  linecolor rgb "black"
 805         set style line 3 linetype 3  linecolor rgb "green"
 806         set style increment user
 807         set ylabel "Indel count"
 808         set xlabel "Read Cycle"
 809         set title "$$args{title}"
 810         plot '-' w l ti 'Insertions', '' w l ti 'Deletions'
 811     ];
 812     for my $len (@{$$opts{dat}{IC}}) { print $fh "$$len[0]\t$$len[1]\n"; } print $fh "end\n";
 813     for my $len (@{$$opts{dat}{IC}}) { print $fh "$$len[0]\t$$len[2]\n"; } print $fh "end\n";
 814     close($fh);
 815     plot($$args{gp});
 816 }
 817
 818
 819
 820
 821
 822
 823
 824 sub has_values
 825 {
 826     my ($opts,@tags) = @_;
 827     for my $tag (@tags)
 828     {
 829         my (@lines) = `cat $$opts{bamcheck} | grep ^$tag | wc -l`;
 830         chomp($lines[0]);
 831         if ( $lines[0]<2 ) { return 0; }
 832     }
 833     return 1;
 834 }
 835
 836 sub plot_mismatches_per_cycle_old
 837 {
 838     my ($opts) = @_;
 839
 840     my $args = get_defaults($opts,"$$opts{prefix}mism-per-cycle.png");
 841     my ($nquals)  = `grep ^MPC $$opts{bamcheck} | awk '\$2==1' | sed 's,\\t,\\n,g' | wc -l`;
 842     my ($ncycles) = `grep ^MPC $$opts{bamcheck} | wc -l`;
 843     chomp($nquals);
 844     chomp($ncycles);
 845     $nquals--;
 846     $ncycles--;
 847     my @gr0_15  = (2..17);
 848     my @gr16_30 = (18..32);
 849     my @gr31_n  = (33..$nquals);
 850     my $gr0_15  = '$'. join('+$',@gr0_15);
 851     my $gr16_30 = '$'. join('+$',@gr16_30);
 852     my $gr31_n  = '$'. join('+$',@gr31_n);
 853
 854     open(my $fh,'>',$$args{gp}) or error("$$args{gp}: $!");
 855     print $fh q[
 856         set terminal png size 600,400 truecolor font "DejaVuSansMono,9"
 857         set output "] . $$args{img} . q["
 858
 859         set key left top
 860         set style data histogram
 861         set style histogram rowstacked
 862
 863         set grid back lc rgb "#aaaaaa"
 864         set ylabel "Number of mismatches"
 865         set xlabel "Read Cycle"
 866         set style fill solid border -1
 867         set title "] . $$args{title} . qq["
 868         set xrange [-1:$ncycles]
 869
 870         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'
 871     ];
 872     close($fh);
 873
 874     plot($$args{gp});
 875 }
 876
 877