plot-vcfstats

#!/usr/bin/env perl
#
# Author: petr.danecek@sanger
#
# Dependencies: 
#	- matplotlib
#		http://matplotlib.sourceforge.net
#	- LaTex/xcolor.sty
#		Download .sty.gz LaTeX class style from http://www.ukern.de/tex/xcolor.html,
#		unpack and install system-wide or place elsewhere and make available by
#		setting the TEXINPUTS environment variable (note the colon)
#			export TEXINPUTS=some/dir:
#		The list of the recognised path can be obtained from `kpsepath tex`
#	
#

use strict;
use warnings;
use Carp;
use Storable qw(dclone);

my $opts = parse_params();
parse_vcfstats($opts);
merge_vcfstats($opts) if @{$$opts{vcfstats}} > 1;

if ( $$opts{make_plots} )
{
	init_plots($opts);
	plot_venn_bars($opts);
	plot_counts_by_AF($opts);
    plot_overlap_by_AF($opts);
	plot_concordance_by_AF($opts);
	plot_concordance_by_sample($opts);
	for my $id (file_ids($opts))
	{
		plot_irc_by_AF($opts,$id);
		plot_irc_by_rlen($opts,$id);
		plot_tstv_by_AF($opts,$id);
		plot_tstv_by_QUAL($opts,$id);
		plot_indel_distribution($opts,$id);
		plot_substitutions($opts,$id);
		plot_per_sample_stats($opts,$id);
		plot_DP($opts,$id);
        plot_hwe($opts,$id);
	}
	plot($opts);
}
create_pdf($opts) unless !$$opts{make_pdf};

exit;

#--------------------------------

sub error
{
    my (@msg) = @_;
    if ( scalar @msg ) { confess @msg; }
    die
        "Usage: plot-vcfstats [OPTIONS] file.chk ...\n",
        "       plot-vcfstats -p outdir/ file.chk ...\n",
        "Options:\n",
        "   -m, --merge                         Merge vcfstats files to STDOUT, skip plotting.\n",
        "   -p, --prefix <path>                 The output files prefix, add a slash to create new directory.\n",
        "   -P, --no-PDF                        Skip the PDF creation step.\n",
        "   -r, --rasterize                     Rasterize PDF images for fast rendering.\n",
        "   -s, --sample-names                  Use sample names for xticks rather than numeric IDs.\n",
        "   -t, --title <string>                Identify files by these titles in plots. Can be given multiple times.\n",
        "   -T, --main-title <string>           Main title for the PDF.\n",
        "   -h, -?, --help                      This help message.\n",
        "\n";
}


sub parse_params
{
    $0 =~ s{^.+/}{};
    my $opts = 
	{ 
        pdf_plots  => 1,
        use_sample_names => 0,
		verbose    => 1,
		make_pdf   => 1,
		make_plots => 1,
		merge      => 0,
		args       => join(' ',$0,@ARGV), 
		img_width  => 11/2.54,
		img_height => 10/2.54,
		id2col     => [ 'orange', 'red', 'darkgreen' ],
		tex => 
		{
			slide3v => { height1 => '7cm', height2 => '7cm',  height3 => '4.5cm' },
			slide3h => { width1  => '15cm', width2 => '10cm', width3 => '8cm' },
		},

        # for file version sanity check
        sections =>
        [
            { 
                id=>'ID',    
                header=>'Definition of sets', 
                exp=>"# ID\t[2]id\t[3]tab-separated file names" 
            },
            { 
                id=>'SN',    
                header=>'SN, Summary numbers', 
                exp=>"# SN  [2]id   [3]key  [4]value" 
            },
            { 
                id=>'TSTV',    
                header=>'# TSTV, transition/transversions:', 
                exp=>"# TSTV\t[2]id\t[3]ts\t[4]tv\t[5]ts/tv\t[6]ts (1st ALT)\t[7]tv (1st ALT)\t[8]ts/tv (1st ALT)" 
            },
            { 
                id=>'SiS',   
                header=>'Sis, Singleton stats', 
                exp=>"# SiS\t[2]id\t[3]allele count\t[4]number of SNPs\t[5]number of transitions\t[6]number of transversions\t[7]number of indels\t[8]repeat-consistent\t[9]repeat-inconsistent\t[10]not applicable" 
            },
            { 
                id=>'AF',    
                header=>'AF, Stats by non-reference allele frequency', 
                exp=>"# AF\t[2]id\t[3]allele frequency\t[4]number of SNPs\t[5]number of transitions\t[6]number of transversions\t[7]number of indels\t[8]repeat-consistent\t[9]repeat-inconsistent\t[10]not applicable" 
            },
            { 
                id=>'IDD',   
                header=>'IDD, InDel distribution', 
                exp=>"# IDD\t[2]id\t[3]length (deletions negative)\t[4]count" 
            },
            { 
                id=>'ST',    
                header=>'ST, Substitution types', 
                exp=>"# ST\t[2]id\t[3]type\t[4]count" 
            },
            { 
                id=>'GCsAF', 
                header=>'GCsAF, Genotype concordance by non-reference allele frequency (SNPs)', 
                exp=>"# GCsAF\t[2]id\t[3]allele frequency\t[4]RR Hom matches\t[5]RA Het matches\t[6]AA Hom matches\t[7]RR Hom mismatches\t[8]RA Het mismatches\t[9]AA Hom mismatches\t[10]dosage r-squared\t[11]number of sites" 
            },
            { 
                id=>'GCiAF', 
                header=>'GCiAF, Genotype concordance by non-reference allele frequency (indels)', 
                exp=>"# GCiAF\t[2]id\t[3]allele frequency\t[4]RR Hom matches\t[5]RA Het matches\t[6]AA Hom matches\t[7]RR Hom mismatches\t[8]RA Het mismatches\t[9]AA Hom mismatches\t[10]dosage r-squared\t[11]number of sites" 
            },
            { 
                id=>'NRDs',
                header=>'Non-Reference Discordance (NRD), SNPs', 
                exp=>"# NRDs\t[2]id\t[3]NRD\t[4]Ref/Ref discordance\t[5]Ref/Alt discordance\t[6]Alt/Alt discordance" 
            },
            { 
                id=>'NRDi',
                header=>'Non-Reference Discordance (NRD), indels', 
                exp=>"# NRDi\t[2]id\t[3]NRD\t[4]Ref/Ref discordance\t[5]Ref/Alt discordance\t[6]Alt/Alt discordance" 
            },
            { 
                id=>'GCsS',  
                header=>'GCsS, Genotype concordance by sample (SNPs)', 
                exp=>"# GCsS\t[2]id\t[3]sample\t[4]non-reference discordance rate\t[5]RR Hom matches\t[6]RA Het matches\t[7]AA Hom matches\t[8]RR Hom mismatches\t[9]RA Het mismatches\t[10]AA Hom mismatches" 
            },
            { 
                id=>'GCiS',  
                header=>'GCiS, Genotype concordance by sample (indels)', 
                exp=>"# GCiS\t[2]id\t[3]sample\t[4]non-reference discordance rate\t[5]RR Hom matches\t[6]RA Het matches\t[7]AA Hom matches\t[8]RR Hom mismatches\t[9]RA Het mismatches\t[10]AA Hom mismatches" 
            },
            { 
                id=>'PSC',   
                header=>'PSC, Per-sample counts', 
                exp=>"# PSC\t[2]id\t[3]sample\t[4]nRefHom\t[5]nNonRefHom\t[6]nHets\t[7]nTransitions\t[8]nTransversions\t[9]nIndels\t[10]average depth\t[11]nSingletons" 
            },
            { 
                id=>'PSI',
                header=>'PSI, Per-sample Indels', 
                exp=>"# PSI\t[2]id\t[3]sample\t[4]in-frame\t[5]out-frame\t[6]not applicable\t[7]out/(in+out) ratio" 
            },
            { 
                id=>'DP',    
                header=>'DP, Depth distribution', 
                exp=>"# DP\t[2]id\t[3]bin\t[4]number of genotypes\t[5]fraction of genotypes (%)" 
            },
            { 
                id=>'FS',    
                header=>'FS, Indel frameshifts', 
                exp=>"# FS\t[2]id\t[3]in-frame\t[4]out-frame\t[5]not applicable\t[6]out/(in+out) ratio\t[7]in-frame (1st ALT)\t[8]out-frame (1st ALT)\t[9]not applicable (1st ALT)\t[10]out/(in+out) ratio (1st ALT)"
            },
            { 
                id=>'ICS',
                header=>'ICS, Indel context summary', 
                exp=>"# ICS\t[2]id\t[3]repeat-consistent\t[4]repeat-inconsistent\t[5]not applicable\t[6]c/(c+i) ratio" 
            },
            { 
                id=>'ICL',
                header=>'ICL, Indel context by length', 
                exp=>"# ICL\t[2]id\t[3]length of repeat element\t[4]repeat-consistent deletions)\t[5]repeat-inconsistent deletions\t[6]consistent insertions\t[7]inconsistent insertions\t[8]c/(c+i) ratio"
            },
            { 
                id=>'QUAL',
                header=>'QUAL, Stats by quality', 
                exp=>"# QUAL\t[2]id\t[3]Quality\t[4]number of SNPs\t[5]number of transitions (1st ALT)\t[6]number of transversions (1st ALT)\t[7]number of indels" 
            },
            { 
                id=>'HWE',
                header=>'HWE', 
                exp=>"# HWE\t[2]id\t[3]1st ALT allele frequency\t[4]Number of observations\t[5]25th percentile\t[6]median\t[7]75th percentile",
            },
        ],
	};
    for my $sec (@{$$opts{sections}}) { $$opts{exp}{$$sec{id}} = $$sec{exp}; $$opts{id2sec}{$$sec{id}} = $sec; }
    while (defined(my $arg=shift(@ARGV)))
    {
        if (                 $arg eq '--no-plots' ) { $$opts{make_plots}=0; next; }
        if ( $arg eq '-P' || $arg eq '--no-PDF' ) { $$opts{make_pdf}=0; next; }
        if ( $arg eq '-r' || $arg eq '--rasterize' ) { $$opts{rasterize}=1; $$opts{pdf_plots} = 0; next; }
        if ( $arg eq '-m' || $arg eq '--merge' ) { $$opts{make_plots}=0; $$opts{make_pdf}=0; $$opts{merge}=1; next; }
        if ( $arg eq '-s' || $arg eq '--sample-names' ) { $$opts{use_sample_names}=1; next; }
        if ( $arg eq '-t' || $arg eq '--title' ) { push @{$$opts{titles}},shift(@ARGV); next; }
        if ( $arg eq '-T' || $arg eq '--main-title' ) { $$opts{main_title} = shift(@ARGV); next; }
        if ( $arg eq '-p' || $arg eq '--prefix' ) { $$opts{prefix}=shift(@ARGV); next; }
        if ( $arg eq '-?' || $arg eq '-h' || $arg eq '--help' ) { error(); }
        if ( -e $arg ) { push @{$$opts{vcfstats}},$arg; next; }
        error("Unknown parameter or non-existent file \"$arg\". Run -h for help.\n");
    }
    if ( !exists($$opts{vcfstats}) ) { error(); }
    if ( !exists($$opts{prefix}) )
    {
        if ( !$$opts{merge} ) { error("Expected -p parameter.\n") }
        $$opts{prefix} = './';
    }
    elsif ( $$opts{merge} ) { error("Only one of -p or -m should be given.\n"); }
    if ( $$opts{merge} && @{$$opts{vcfstats}} < 2 ) { error("Nothing to merge\n") }

	$$opts{dir} = $$opts{prefix};
	$$opts{dir} =~ s{/[^/]+$}{/};
	if ( !($$opts{dir}=~m{/$}) ) { $$opts{dir} = './'; }
	if ( !($$opts{prefix}=~m{/$}) && !($$opts{prefix}=~/-$/) ) { $$opts{prefix} .= '-'; }
	$$opts{lprefix} = $$opts{prefix};
	$$opts{lprefix} =~ s{^.+/}{};
	$$opts{logfile} = "$$opts{lprefix}plot-vcfstats.log";
	if ( !-d $$opts{dir} ) { `mkdir -p $$opts{dir}`; } 
	`> $$opts{dir}$$opts{logfile}` unless $$opts{merge};
    return $opts;
}


sub plot
{
    my ($opts) = @_;
	if ( !exists($$opts{plt_fh}) ) { return; }
	close($$opts{plt_fh}) or error("close $$opts{plt_fh}");
    my $cmd = "python $$opts{plt_file}";
	print STDERR "Plotting graphs: $cmd\n" unless !$$opts{verbose};
    system($cmd);
    if ( $? ) { error("The command exited with non-zero status $?:\n\t$cmd\n\n"); }
}


sub parse_vcfstats
{
    my ($opts) = @_;
    for (my $i=0; $i<@{$$opts{vcfstats}}; $i++) { parse_vcfstats1($opts,$i); }

    # Check sanity
    if ( !exists($$opts{dat}{ID}{0}) )
    {
        error("Sanity check failed: no stats found by vcfstats??\n");
    }

	# Set titles
	my %file2title;
    my %title2file;
	if ( exists($$opts{titles}) )
	{
		for (my $i=0; $i<@{$$opts{titles}}; $i++)
		{
			if ( !exists($$opts{dat}{ID}{$i}) ) { next; }
			$file2title{$$opts{dat}{ID}{$i}[0][0]} = $$opts{titles}[$i];
            $title2file{$$opts{titles}[$i]} = $$opts{dat}{ID}{$i}[0][0];
		}
	}
	for my $id (file_ids($opts))
	{
        if ( @{$$opts{dat}{ID}{$id}[0]}>1 ) { next; }   # shared stats of two files
        my $file = $$opts{dat}{ID}{$id}[0][0];
        if ( !exists($file2title{$file}) )  # create short title
        {
            my $bname = $file;
            $bname =~ s{^.*/}{};
            $bname =~ s{\.vcf\.gz$}{}i;
            if ( length($bname) > 5 ) { $bname = substr($bname,0,5); }
            my $i = 0;
            my $title = $bname;
            while ( exists($title2file{$title}) ) { $title = $bname.chr(66+$i); $i++; }
            $file2title{$file} = $title;
            $title2file{$title} = $file;
        }
	}
	for my $id (file_ids($opts))
	{
		my @titles;
		for my $file (@{$$opts{dat}{ID}{$id}[0]}) { push @titles, $file2title{$file} if $file2title{$file}; }
        $$opts{title}{$id} = join(' + ',@titles);
	}

	# mapping from file names to list of IDs
	for my $id (file_ids($opts))
	{
		my @files;
		for my $file (@{$$opts{dat}{ID}{$id}[0]})
		{
			push @{$$opts{file2ids}{$file}}, $id;
		}
	}

    # check sanity of the file merge: were the correct files merged?
    if ( exists($$opts{coalesced_files}) && $$opts{verbose} )
    {
        print STDERR "The vcfstats outputs have been merged as follows:\n";
        my %printed;
        for my $id (keys %{$$opts{coalesced_files}})
        {
            for (my $i=0; $i<@{$$opts{coalesced_files}{$id}}; $i++)
            {
                for (my $j=0; $j<@{$$opts{coalesced_files}{$id}[$i]}; $j++)
                {
                    if ( exists($printed{$$opts{dat}{ID}{$id}[$i][$j]}) ) { next; }
                    print STDERR "\t$$opts{dat}{ID}{$id}[$i][$j]\n";
                    for my $file (keys %{$$opts{coalesced_files}{$id}[$i][$j]})
                    {
                        my $n = $$opts{coalesced_files}{$id}[$i][$j]{$file};
                        print STDERR "\t\t$file", ($n>1 ? "\t..\t${n}x" : ''),"\n";
                    }
                    $printed{$$opts{dat}{ID}{$id}[$i][$j]} = 1;
                }
            }
        }
    }
}

sub add_to_values
{
    my ($dst,$src,$cmp) = @_;
    my $id = 0;
    my $is = 0;
    while ($is<@$src)
    {
        while ( $id<@$dst && &$cmp($$src[$is][0],$$dst[$id][0])>0 ) { $id++; }
        if ( $id<@$dst && !&$cmp($$src[$is][0],$$dst[$id][0]) ) 
        { 
            for (my $j=1; $j<@{$$src[$is]}; $j++) { $$dst[$id][$j] += $$src[$is][$j]; }
        }
        else { splice(@$dst,$id,0,$$src[$is]); }
        $is++;
    }
}

sub add_to_sample_values
{
    my ($dst,$src) = @_;
    my %id2i;
    for (my $i=0; $i<@$dst; $i++)
    {
        $id2i{$$dst[$i][0]} = $i;
    }
    for (my $i=0; $i<@$src; $i++)
    {
        if ( !exists($id2i{$$src[$i][0]}) ) { error("Whoops, no such dst sample: $$src[$i][0]\n"); }
        my $di = $id2i{$$src[$i][0]};
        for (my $j=1; $j<@{$$src[$i]}; $j++)
        {
            $$dst[$di][$j] += $$src[$i][$j];
        }
    }
}

sub merge_PSC
{
    my ($a,$b,$n) = @_;
    for (my $i=0; $i<@$a; $i++) { $$a[$i][7] *= $n; }   # average DP
    add_to_sample_values($a,$b);
    for (my $i=0; $i<@$a; $i++) { $$a[$i][7] /= $n+1; }
}

sub merge_PSI
{
    my ($a,$b,$n) = @_;
    add_to_sample_values($a,$b);
    for (my $i=0; $i<@$b; $i++) { $$a[$i][4] = sprintf "%.2f", ($$a[$i][1]+$$a[$i][2] ? $$a[$i][2]/($$a[$i][1]+$$a[$i][2]) : 0); }
}

sub rglob
{
    my ($a,$b) = @_;
    if ( $a eq $b ) { return $a; }
    $a =~ s/\*//;
    my $la = length($a);
    my $lb = length($b);
    my $i = 0;
    while ( $i<$la && $i<$lb && substr($a,$i,1) eq substr($b,$i,1) ) { $i++; }
    $la--; $lb--;
    while ( $la>$i && $lb>$i && substr($a,$la,1) eq substr($b,$lb,1) ) { $la--; $lb--; }
    $la = $la==$i && $lb==$i ? 1 : $la-$i;
    substr($a,$i,$la,'*');
    return $a;
}

sub merge_id
{
    # merge id filenames
    my ($opts,$dst,$src,$id) = @_;
    for (my $i=0; $i<@{$$src{$id}}; $i++)
    {
        for (my $j=0; $j<@{$$src{$id}[$i]}; $j++)
        {
            my $gname = rglob($$dst{$id}[$i][$j], $$src{$id}[$i][$j]);
            $$dst{$id}[$i][$j] = $gname;
            $$opts{coalesced_files}{$id}[$i][$j]{$$src{$id}[$i][$j]}++;
        }
    }
}

sub add_to_avg
{
    my ($dst,$src,$n) = @_;
    for (my $i=0; $i<@$src; $i++)
    {
        if ( ref($$dst[$i]) eq 'ARRAY' )
        {
            for (my $j=0; $j<@{$$dst[$i]}; $j++)
            {
                $$dst[$i][$j] = ($n*$$dst[$i][$j]+$$src[$i][$j])/($n+1);
            }
        }
        else
        {
            $$dst[$i] = ($n*$$dst[$i]+$$src[$i])/($n+1);
        }
    }
}

sub cmp_str($$) { my ($a,$b) = @_; return $a cmp $b; }
sub cmp_num($$) { my ($a,$b) = @_; return $a <=> $b; }
sub cmp_num_op($$) 
{ 
    # numeric compare with operators
    # Cases like <3, >500 make it complicated
    my ($a,$b) = @_;
    my $xa = '=';
    my $xb = '=';
    if ( $a=~/^(\D+)/ ) { $xa = $1; $a = $'; }
    if ( $b=~/^(\D+)/ ) { $xb = $1; $b = $'; }
    if ( $a==$b ) { return $xa cmp $xb; }
    $a <=> $b;
}

sub merge_dp
{
    my ($a,$b) = @_;
    add_to_values($a,$b,\&cmp_num_op);
    # recalculate fraction of GTs, cannot be simply summed
    my $sum = 0;
    for (my $i=0; $i<@$a; $i++) { $sum += $$a[$i][1]; }
    for (my $i=0; $i<@$a; $i++) { $$a[$i][2] = $$a[$i][1]*100./$sum; }
}

sub merge_GCsS
{
    my ($a,$b,$n) = @_;
    # average the non-ref discordance rate
    for (my $i=0; $i<@$a; $i++) { $$a[$i][1] *= $n; }
    add_to_sample_values($a,$b);
    for (my $i=0; $i<@$a; $i++) { $$a[$i][1] /= $n+1; }
}

sub merge_FS
{
    my ($a,$b) = @_;
    for (my $i=0; $i<@$a; $i++)
    {
        for (my $j=0; $j<3; $j++) { $$a[$i][$j] += $$b[$i][$j]; }
        $$a[$i][3] = sprintf "%.2f", ($$a[$i][0] + $$a[$i][1]) ? $$a[$i][1]/($$a[$i][0] + $$a[$i][1]) : 0;

        for (my $j=4; $j<7; $j++) { $$a[$i][$j] += $$b[$i][$j]; }
        $$a[$i][7] = sprintf "%.2f", ($$a[$i][4] + $$a[$i][5]) ? $$a[$i][5]/($$a[$i][4] + $$a[$i][5]) : 0;
    }
}

sub merge_ICS
{
    my ($a,$b) = @_;
    for (my $i=0; $i<@$a; $i++)
    {
        for (my $j=0; $j<3; $j++) { $$a[$i][$j] += $$b[$i][$j]; }
        $$a[$i][3] = sprintf "%.4f", ($$a[$i][0] + $$a[$i][1]) ? $$a[$i][0]/($$a[$i][0] + $$a[$i][1]) : 0;
    }
}

sub merge_ICL
{
    my ($a,$b) = @_;
    for (my $i=0; $i<@$a; $i++)
    {
        for (my $j=1; $j<5; $j++) { $$a[$i][$j] += $$b[$i][$j]; }
        $$a[$i][5] = sprintf "%.4f", ($$a[$i][2] + $$a[$i][4]) ? ($$a[$i][1] + $$a[$i][3])/($$a[$i][1] + $$a[$i][2] + $$a[$i][3] + $$a[$i][4]) : 0;
    }
}

sub merge_TSTV
{
    my ($a,$b) = @_;
    for (my $i=0; $i<@$a; $i++)
    {
        for (my $j=0; $j<2; $j++) { $$a[$i][$j] += $$b[$i][$j]; }
        $$a[$i][2] = sprintf "%.2f", $$a[$i][1] ? $$a[$i][0]/$$a[$i][1] : 0;

        for (my $j=3; $j<5; $j++) { $$a[$i][$j] += $$b[$i][$j]; }
        $$a[$i][5] = sprintf "%.2f", $$a[$i][4] ? $$a[$i][3]/$$a[$i][4] : 0;
    }
}

sub merge_GCsAF
{
    my ($a,$b) = @_;
    # recalculate r2
    for (my $i=0; $i<@$a; $i++) { $$a[$i][7] *= $$a[$i][8]; }
    for (my $i=0; $i<@$b; $i++) { $$b[$i][7] *= $$b[$i][8]; }
    add_to_values($a,$b,\&cmp_num_op);
    for (my $i=0; $i<@$a; $i++) { $$a[$i][7] /= $$a[$i][8]; }
}

sub parse_vcfstats1
{
    my ($opts,$i) = @_;
    my $file = $$opts{vcfstats}[$i];
	print STDERR "Parsing vcfstats output: $file\n" unless !$$opts{verbose};
    open(my $fh,'<',$file) or error("$file: $!");
    my $line = <$fh>;
    if ( !$line or !($line=~/^# This file was produced by \S*vcfstats/) ) { error("Sanity check failed: was this file generated by vcfstats?"); }
    my %dat;
    while ($line=<$fh>)
    {
        $line =~ s/\s*$//;
        if ( $line=~/^#\s+(\S+)\t/ ) 
        {
            $$opts{def_line}{$1} = $line;
            next; 
        }
        if ( $line=~/^#/ ) { next; }
        my @items = split(/\t/,$line);
        if ( $items[0] eq 'SN' )
        {
            $dat{$items[1]}{$items[2]} = splice(@items,3);
            next;
        }
        push @{$dat{$items[0]}{$items[1]}}, [splice(@items,2)];
    }
    close($fh);
    for my $a (keys %dat)
    {
        if ( !exists($$opts{dat}{$a}) ) { $$opts{dat}{$a} = $dat{$a}; next; } # first vcfstats file
        for my $b (keys %{$dat{$a}})
        {
            # Merging multiple vcfstats files. Honestly, this is quite hacky.
            if ( !exists($$opts{dat}{$a}{$b}) ) { $$opts{dat}{$a}{$b} = $dat{$a}{$b}; next; } # copy all, first occurance

            if ( $a eq 'ID' ) { merge_id($opts,$$opts{dat}{$a},$dat{$a},$b); }
            elsif ( ref($dat{$a}{$b}) ne 'ARRAY' ) {$$opts{dat}{$a}{$b} += $dat{$a}{$b} unless $b eq 'number of samples:'; } # SN, Summary numbers, do not sum sample counts
            elsif ( $a eq 'NRDs' ) { add_to_avg($$opts{dat}{$a}{$b},$dat{$a}{$b},$i); }
            elsif ( $a eq 'NRDi' ) { add_to_avg($$opts{dat}{$a}{$b},$dat{$a}{$b},$i); }
            elsif ( $a eq 'DP' ) { merge_dp($$opts{dat}{$a}{$b},$dat{$a}{$b}); }
            elsif ( $a eq 'GCsS' ) { merge_GCsS($$opts{dat}{$a}{$b},$dat{$a}{$b},$i); }
            elsif ( $a eq 'GCiS' ) { merge_GCsS($$opts{dat}{$a}{$b},$dat{$a}{$b},$i); }
            elsif ( $a eq 'GCsAF' ) { merge_GCsAF($$opts{dat}{$a}{$b},$dat{$a}{$b},$i); }
            elsif ( $a eq 'GCiAF' ) { merge_GCsAF($$opts{dat}{$a}{$b},$dat{$a}{$b},$i); }
            elsif ( $a eq 'ST' ) { add_to_values($$opts{dat}{$a}{$b},$dat{$a}{$b},\&cmp_str); }
            elsif ( $a eq 'PSC') { merge_PSC($$opts{dat}{$a}{$b},$dat{$a}{$b},$i); }
            elsif ( $a eq 'PSI') { merge_PSI($$opts{dat}{$a}{$b},$dat{$a}{$b},$i); }
            elsif ( $a eq 'IDD') { add_to_values($$opts{dat}{$a}{$b},$dat{$a}{$b},\&cmp_num); }
            elsif ( $a eq 'FS') { merge_FS($$opts{dat}{$a}{$b},$dat{$a}{$b}); }
            elsif ( $a eq 'ICS') { merge_ICS($$opts{dat}{$a}{$b},$dat{$a}{$b}); }
            elsif ( $a eq 'ICL') { merge_ICL($$opts{dat}{$a}{$b},$dat{$a}{$b}); }
            elsif ( $a eq 'TSTV') { merge_TSTV($$opts{dat}{$a}{$b},$dat{$a}{$b},$i); }
            elsif ( $a eq 'DBG' ) { next; }
            else { add_to_values($$opts{dat}{$a}{$b},$dat{$a}{$b},\&cmp_num_op); } # SiS AF IDD
        }
    }
}

sub file_ids
{
	my ($opts) = @_;
	my $id = 0;
	my @out;
	while ( exists($$opts{dat}{ID}) && exists($$opts{dat}{ID}{$id}) ) { push @out, $id++; }
	return @out;
}

sub tprint
{
	my ($fh,@txt) = @_;
	for my $txt (@txt)
	{
		$txt =~ s/\n[ \t]+/\n/g;        # eat leading tabs
        while ( ($txt =~ /\n\t*\\t/) )
        {
            $txt =~ s/(\n\t*)\\t/$1\t/g;    # replace escaped tabs (\\t) with tabs
        }
		print $fh $txt;
	}
}

sub init_plots
{
    my ($opts) = @_;

	$$opts{plt_file} = "$$opts{prefix}plot.py";

    my $titles = "# Title abbreviations:\n";
    for my $id (file_ids($opts))
    {
        $titles .= "# \t $id .. $$opts{title}{$id} .. $$opts{dat}{ID}{$id}[0][0]\n";
    }
    $titles .= "#";

	open(my $fh,'>',$$opts{plt_file}) or error("$$opts{plt_file}: $!");
	tprint $fh, qq[
		# This file was produced by plot-vcfstats, the command line was:
		#	$$opts{args}
		#
		# Edit as necessary and recreate the plots by running
		#	python $$opts{plt_file} 
		#
        $titles

        # Set to 1 to plot in PDF instead of PNG
        pdf_plots = $$opts{pdf_plots}

        # Set to 1 to use sample names for xticks instead of numeric sequential IDs
        #   and adjust margins and font properties if necessary
        sample_names   = $$opts{use_sample_names}
        sample_margins = {'right':0.98, 'left':0.07, 'bottom':0.2}
        sample_font    = {'rotation':45, 'ha':'right', 'fontsize':8}

        if sample_names==0: sample_margins=(); sample_font=();

		import matplotlib as mpl
		mpl.use('Agg')
		import matplotlib.pyplot as plt

        import csv
        csv.register_dialect('tab', delimiter='\\t', quoting=csv.QUOTE_NONE)

        import numpy
        def smooth(x,window_len=11,window='hanning'):
        \\tif x.ndim != 1: raise ValueError, "smooth only accepts 1 dimension arrays."
        \\tif x.size < window_len: return x
        \\tif window_len<3: return x
        \\tif not window in ['flat', 'hanning', 'hamming', 'bartlett', 'blackman']: raise ValueError, "Window is on of 'flat', 'hanning', 'hamming', 'bartlett', 'blackman'"
        \\ts = numpy.r_[x[window_len-1:0:-1],x,x[-1:-window_len:-1]]
        \\tif window == 'flat': # moving average
        \\t\\tw = numpy.ones(window_len,'d')
        \\telse:
        \\t\\tw = eval('numpy.'+window+'(window_len)')
        \\ty = numpy.convolve(w/w.sum(),s,mode='valid')
        \\treturn y[(window_len/2-1):-(window_len/2)]

		];
	$$opts{plt_fh} = $fh;
}

sub percentile
{
    my ($p,@vals) = @_;
    my $N = 0;
    for my $val (@vals) { $N += $val; }
    my $n = $p*($N+1)/100.;
    my $k = int($n);
    my $d = $n-$k;
    if ( $k<=0 ) { return 0; }
    if ( $k>=$N ) { return scalar @vals-1; }
    my $cnt;
    for (my $i=0; $i<@vals; $i++)
    { 
        $cnt += $vals[$i]; 
        if ( $cnt>=$k ) { return $i; }
    }
    error("FIXME: this should not happen [percentile]\n");
}

sub get_values
{
	my ($opts,$id,$key,$i,$j) = @_;
	if ( !exists($$opts{dat}{$key}) ) { return (); }
	if ( !exists($$opts{dat}{$key}{$id}) ) { return (); }
    my $fields_ok = 1;
    if ( !exists($$opts{exp}{$key}) ) { error("todo: sanity check for $key\n"); } 
    if ( exists($$opts{def_line}{$key}) && $$opts{def_line}{$key} ne $$opts{exp}{$key} && !$$opts{def_line_warned}{$key} )
    {
        warn("Warning: Possible version mismatch, the definition line differs\n\texpected: $$opts{exp}{$key}\n\tfound:    $$opts{def_line}{$key}\n");
        $$opts{def_line_warned}{$key} = 1;
    }
    if ( defined $i )
    {
        if ( defined $j ) { return $$opts{dat}{$key}{$id}[$i][$j]; }
        return (@{$$opts{dat}{$key}{$id}[$i]});
    }
	return (@{$$opts{dat}{$key}{$id}});
}

sub get_value
{
	my ($opts,$id,$key) = @_;
	if ( !exists($$opts{dat}{$id}) ) { return undef; }
	if ( !exists($$opts{dat}{$id}{$key}) ) { return undef}
	return $$opts{dat}{$id}{$key};
}

sub plot_venn_bars
{
    my ($opts) = @_;

	my @ids  = file_ids($opts);
	if ( @ids != 3 ) { return; }

	my (@snps,@indels,@tstv,@snp_titles,@indel_titles);
	for my $id (0..2)
	{
		push @snps, get_value($opts,$id,'number of SNPs:');
		push @indels, get_value($opts,$id,'number of indels:');
		push @tstv, sprintf("%.2f",get_values($opts,$id,'TSTV',0,5));
		push @snp_titles, "$$opts{title}{$id}\\nts/tv $tstv[$id]\\n" .bignum($snps[$id]);
		my @fs = get_values($opts,$id,'FS');
        my $fs = @fs ? "frm $fs[0][3]\\n" : '';
		push @indel_titles, "$$opts{title}{$id}\\n$fs" .bignum($indels[$id]);
	}

    my $fh = $$opts{plt_fh};
	tprint $fh, qq[
			fig = plt.figure(figsize=($$opts{img_width},$$opts{img_height}))
	 		ax1 = fig.add_subplot(111)
	 		ax1.bar([1,2,3],[$snps[0],$snps[2],$snps[1]],align='center',color='$$opts{id2col}[0]',width=0.3)
			ax1.ticklabel_format(style='sci', scilimits=(0,0), axis='y') 
			ax1.set_xlim(0.5,3.5)
			plt.xticks([1,2,3],('$snp_titles[0]','$snp_titles[2]','$snp_titles[1]'))
			plt.title('Number of SNPs')
			plt.subplots_adjust(right=0.95,bottom=0.15)
	 		plt.savefig('$$opts{prefix}venn_bars.snps.png')
	 		if pdf_plots: plt.savefig('$$opts{prefix}venn_bars.snps.pdf')
	 		plt.close()

			fig = plt.figure(figsize=($$opts{img_width},$$opts{img_height}))
	 		ax1 = fig.add_subplot(111)
	 		ax1.bar([1,2,3],[$indels[0],$indels[2],$indels[1]],align='center',color='$$opts{id2col}[1]',width=0.3)
			ax1.ticklabel_format(style='sci', scilimits=(0,0), axis='y') 
			ax1.set_xlim(0.5,3.5)
			plt.xticks([1,2,3],('$indel_titles[0]','$indel_titles[2]','$indel_titles[1]'))
			plt.title('Number of indels')
			plt.subplots_adjust(right=0.95,bottom=0.15)
	 		plt.savefig('$$opts{prefix}venn_bars.indels.png')
	 		if pdf_plots: plt.savefig('$$opts{prefix}venn_bars.indels.pdf')
	 		plt.close()
		];
}

sub plot_per_sample_stats
{
	my ($opts,$id) = @_;
	my @vals = get_values($opts,$id,'PSC');
	if ( !@vals ) { return; }

	my $fh   = $$opts{plt_fh};
	my $img  = "$$opts{prefix}tstv_by_sample.$id";
	my $img2 = "$$opts{prefix}hets_by_sample.$id";
	my $img3 = "$$opts{prefix}snps_by_sample.$id";
	my $img4 = "$$opts{prefix}indels_by_sample.$id";
	my $img5 = "$$opts{prefix}singletons_by_sample.$id";
	my $img6 = "$$opts{prefix}dp_by_sample.$id";

	tprint $fh, "
			dat = [
		";
	for (my $i=0; $i<@vals; $i++)
	{ 
		my $tstv = $vals[$i][5] ? $vals[$i][4]/$vals[$i][5] : 0;
        my $hethom = $vals[$i][2] ? $vals[$i][3]/$vals[$i][2] : 0;
		tprint $fh, sprintf("\t[ %d, %f, %f, %d, %d, %f, %d, r'%s' ],\n", $i, $tstv, $hethom, $vals[$i][4]+$vals[$i][5], $vals[$i][6], $vals[$i][7], $vals[$i][8], $vals[$i][0]);
	}
	tprint $fh, "]

			fig = plt.figure(figsize=(2*$$opts{img_width},$$opts{img_height}*0.7))
			ax1 = fig.add_subplot(111)
			ax1.plot([row[0] for row in dat], [row[1] for row in dat], 'o', color='$$opts{id2col}[$id]',mec='$$opts{id2col}[$id]')
			ax1.set_ylabel('Ts/Tv')
			ax1.set_ylim(min(row[1] for row in dat)-0.1,max(row[1] for row in dat)+0.1)
            if sample_names: 
            \\t     plt.xticks([row[0] for row in dat],[row[7] for row in dat],**sample_font)
			\\t     plt.subplots_adjust(**sample_margins)
            else:
			\\t     plt.subplots_adjust(right=0.98,left=0.07,bottom=0.17)
			\\t     ax1.set_xlabel('Sample ID')
			plt.title('$$opts{title}{$id}')
			plt.savefig('$img.png')
			if pdf_plots: plt.savefig('$img.pdf')
			plt.close()

			fig = plt.figure(figsize=(2*$$opts{img_width},$$opts{img_height}*0.7))
			ax1 = fig.add_subplot(111)
			ax1.plot([row[0] for row in dat], [row[2] for row in dat], 'o', color='$$opts{id2col}[$id]',mec='$$opts{id2col}[$id]')
			ax1.set_ylabel('nHet(RA) / nHom(AA)')
			ax1.ticklabel_format(style='sci', scilimits=(0,0), axis='y')
            if sample_names: 
            \\t     plt.xticks([row[0] for row in dat],[row[7] for row in dat],**sample_font)
			\\t     plt.subplots_adjust(**sample_margins)
            else:
			\\t     plt.subplots_adjust(right=0.98,left=0.07,bottom=0.17)
			\\t     ax1.set_xlabel('Sample ID')
			plt.title('$$opts{title}{$id}')
			plt.savefig('$img2.png')
			if pdf_plots: plt.savefig('$img2.pdf')
			plt.close()

			fig = plt.figure(figsize=(2*$$opts{img_width},$$opts{img_height}*0.7))
			ax1 = fig.add_subplot(111)
			ax1.plot([row[0] for row in dat], [row[3] for row in dat], 'o', color='$$opts{id2col}[$id]',mec='$$opts{id2col}[$id]')
			ax1.set_ylabel('Number of SNPs')
			ax1.ticklabel_format(style='sci', scilimits=(0,0), axis='y')
            if sample_names: 
            \\t     plt.xticks([row[0] for row in dat],[row[7] for row in dat],**sample_font)
			\\t     plt.subplots_adjust(**sample_margins)
            else:
			\\t     plt.subplots_adjust(right=0.98,left=0.07,bottom=0.17)
			\\t     ax1.set_xlabel('Sample ID')
			plt.title('$$opts{title}{$id}')
			plt.savefig('$img3.png')
			if pdf_plots: plt.savefig('$img3.pdf')
			plt.close()

			fig = plt.figure(figsize=(2*$$opts{img_width},$$opts{img_height}*0.7))
			ax1 = fig.add_subplot(111)
			ax1.plot([row[0] for row in dat], [row[4] for row in dat], 'o', color='$$opts{id2col}[$id]',mec='$$opts{id2col}[$id]')
			ax1.set_ylabel('Number of indels')
			ax1.ticklabel_format(style='sci', scilimits=(0,0), axis='y')
            if sample_names: 
            \\t     plt.xticks([row[0] for row in dat],[row[7] for row in dat],**sample_font)
			\\t     plt.subplots_adjust(**sample_margins)
            else:
			\\t     plt.subplots_adjust(right=0.98,left=0.07,bottom=0.17)
			\\t     ax1.set_xlabel('Sample ID')
			plt.title('$$opts{title}{$id}')
			plt.savefig('$img4.png')
			if pdf_plots: plt.savefig('$img4.pdf')
			plt.close()

			fig = plt.figure(figsize=(2*$$opts{img_width},$$opts{img_height}*0.7))
			ax1 = fig.add_subplot(111)
			ax1.plot([row[0] for row in dat], [row[6] for row in dat], 'o', color='$$opts{id2col}[$id]',mec='$$opts{id2col}[$id]')
			ax1.set_ylabel('Number of singletons')
			ax1.ticklabel_format(style='sci', scilimits=(0,0), axis='y')
            if sample_names: 
            \\t     plt.xticks([row[0] for row in dat],[row[7] for row in dat],**sample_font)
			\\t     plt.subplots_adjust(**sample_margins)
            else:
			\\t     plt.subplots_adjust(right=0.98,left=0.07,bottom=0.17)
			\\t     ax1.set_xlabel('Sample ID')
			plt.title('$$opts{title}{$id}')
			plt.savefig('$img5.png')
			if pdf_plots: plt.savefig('$img5.pdf')
			plt.close()

			fig = plt.figure(figsize=(2*$$opts{img_width},$$opts{img_height}*0.7))
			ax1 = fig.add_subplot(111)
			ax1.plot([row[0] for row in dat], [row[5] for row in dat], 'o', color='$$opts{id2col}[$id]',mec='$$opts{id2col}[$id]')
			ax1.set_ylabel('Average depth')
			ax1.ticklabel_format(style='sci', scilimits=(0,0), axis='y')
            if sample_names: 
            \\t     plt.xticks([row[0] for row in dat],[row[7] for row in dat],**sample_font)
			\\t     plt.subplots_adjust(**sample_margins)
            else:
			\\t     plt.subplots_adjust(right=0.98,left=0.07,bottom=0.17)
			\\t     ax1.set_xlabel('Sample ID')
			plt.title('$$opts{title}{$id}')
			plt.savefig('$img6.png')
			if pdf_plots: plt.savefig('$img6.pdf')
			plt.close()
		";
}

sub plot_DP
{
	my ($opts,$id) = @_;
	my @vals = get_values($opts,$id,'DP');
	if ( !@vals ) { return; }

	my $fh   = $$opts{plt_fh};
	my $img  = "$$opts{prefix}depth.$id";

	tprint $fh, "
			dat = [
		";
    my $sum = 0;
	for (my $i=0; $i<@vals; $i++)
	{ 
        if ( $sum>99. ) { last; }
        if ( !($vals[$i][0]=~/^\d+$/) ) { next; }  # DP ">500" case
        $sum += $vals[$i][2];
		tprint $fh, sprintf("\t[ %d, %f, %f ],\n", $vals[$i][0], $sum, $vals[$i][2]);
	}
	tprint $fh, "]

			fig = plt.figure(figsize=($$opts{img_width},$$opts{img_height}))
			ax1 = fig.add_subplot(111)
            ax1.plot([row[0] for row in dat], [row[2] for row in dat], '-^', color='k')
            ax1.set_ylabel('Number of genotypes [%]',color='k')
			ax1.set_xlabel('Depth')
            ax2 = ax1.twinx()
			ax2.plot([row[0] for row in dat], [row[1] for row in dat], '-o', color='$$opts{id2col}[$id]')
			ax2.set_ylabel('Cumulative number of genotypes [%]',color='$$opts{id2col}[$id]')
            for tl in ax2.get_yticklabels(): tl.set_color('$$opts{id2col}[$id]')
			plt.subplots_adjust(left=0.15,bottom=0.15,right=0.87)
			plt.title('$$opts{title}{$id}')
			plt.savefig('$img.png')
			if pdf_plots: plt.savefig('$img.pdf')
			plt.close()
		";
}

sub plot_hwe
{
	my ($opts,$id) = @_;
	my @vals = get_values($opts,$id,'HWE');
	if ( !@vals ) { return; }

	my $fh   = $$opts{plt_fh};
	my $img  = "$$opts{prefix}hwe.$id";

	tprint $fh, "
			dat = [
		";
    my $sum = 0;
	for (my $i=0; $i<@vals; $i++)
	{ 
        if ( !$vals[$i][1] ) { next; }
		tprint $fh, "\t[", join(',', @{$vals[$i]}), "],\n";
	}
	tprint $fh, "]
            x  = [row[0] for row in dat]
            y1 = smooth(numpy.array([row[2] for row in dat]),40,'hanning')
            y2 = smooth(numpy.array([row[3] for row in dat]),40,'hanning')
            y3 = smooth(numpy.array([row[4] for row in dat]),40,'hanning')
            dp = smooth(numpy.array([row[1] for row in dat]),40,'hanning')
            hwe = []
            for af in x: hwe.append(2*af/100.*(1-af/100.))

			fig = plt.figure(figsize=($$opts{img_width},$$opts{img_height}))
            ax1 = fig.add_subplot(111)
            plots  = ax1.plot(x,hwe,'--',color='#ff9900',label='HWE')
            plots += ax1.plot(x,y2,color='#ff9900',label='Median')
            plots += ax1.plot(x,y3,color='#ffe0b2',label='25-75th percentile')
            ax1.fill_between(x,y1,y3, facecolor='#ffeacc',edgecolor='#ffe0b2')
            ax1.set_ylabel('Fraction of hets',color='#ff9900')
            ax1.set_xlabel('Allele frequency [%]')
            for tl in ax1.get_yticklabels(): tl.set_color('#ff9900')
            ax2 = ax1.twinx()
            plots += ax2.plot(x,dp, 'k', label='Number of sites')
            ax2.set_ylabel('Number of sites')
            ax2.set_yscale('log')
            labels = [l.get_label() for l in plots]
            plt.legend(plots,labels,numpoints=1,markerscale=2,loc='center',prop={'size':9},frameon=False)
            plt.subplots_adjust(left=0.15,bottom=0.15,right=0.86)
			plt.title('$$opts{title}{$id}')
			plt.savefig('$img.png')
			if pdf_plots: plt.savefig('$img.pdf')
			plt.close()
		";
}

sub plot_irc_by_rlen
{
    my ($opts,$id) = @_;
	my @vals = get_values($opts,$id,'ICL');
    if ( !@vals ) { return; }

    my $fh  = $$opts{plt_fh};
	my $img = "$$opts{prefix}irc_by_rlen.$id";

    my $has_vals = 0;
	tprint $fh, "
			dat = [
		";
	for (my $i=0; $i<@vals; $i++)
	{ 
		tprint $fh, sprintf("\t[ %d, %d, %f, %f ],\n", 
            $vals[$i][0],                                                   # repeat type
            $vals[$i][1]+$vals[$i][2]+$vals[$i][3]+$vals[$i][4],            # total count
            $vals[$i][2] ? $vals[$i][1]/($vals[$i][1]+$vals[$i][2]) : 0,    # IRC of deletions
            $vals[$i][4] ? $vals[$i][3]/($vals[$i][3]+$vals[$i][4]) : 0     # IRC of insertions
            );
	}
	tprint $fh, "]

			fig = plt.figure(figsize=($$opts{img_width},$$opts{img_height}))
			ax1 = fig.add_subplot(111)
            ax2 = ax1.twinx()
			plots =  ax1.plot([row[0] for row in dat], [row[2] for row in dat], 'o-', ms=3, color='r', mec='r', label='Deletions'); 
			plots += ax1.plot([row[0] for row in dat], [row[3] for row in dat], 'o-', ms=3, color='g', mec='g', label='Insertions');
			plots += ax2.plot([row[0] for row in dat], [row[1] for row in dat], '^--', ms=3, color='k', mec='k', label='Number of sites');

            ax2.ticklabel_format(style='sci',scilimits=(-3,4))
            ax1.set_xlabel('Repeat type');
            ax1.set_ylabel('Repeat Consistency');
            ax2.set_ylabel('Number of sites');
            labels = [l.get_label() for l in plots]
            ax1.legend(plots,labels,numpoints=1,markerscale=1,loc='best',prop={'size':9},frameon=False)

			plt.subplots_adjust(right=0.88,left=0.15,bottom=0.11)
			plt.title('$$opts{title}{$id}')
			plt.savefig('$img.png')
			if pdf_plots: plt.savefig('$img.pdf')
			plt.close()
		";
}

sub plot_irc_by_AF
{
    my ($opts,$id) = @_;
	my @vals = get_values($opts,$id,'AF');
    if ( !@vals ) { return; }

    my $fh   = $$opts{plt_fh};
	my $img  = "$$opts{prefix}irc_by_af.$id";
	my $vals = rebin_values(\@vals,8,0);

    my $has_vals = 0;
	tprint $fh, "
			dat = [
		";
	for (my $i=0; $i<@$vals; $i++)
	{ 
		if ( $$vals[$i][5] + $$vals[$i][6] == 0 ) { next; }
        $has_vals++;
		tprint $fh, sprintf("\t[ %f, %d, %f, %f ],\n", 
			$$vals[$i][0]/100., 
			$$vals[$i][5] + $$vals[$i][6],
			$$vals[$i][6] ? $$vals[$i][5]/($$vals[$i][5] + $$vals[$i][6]) : 0,
			($$vals[$i][6] + $$vals[$i][7]) ? $$vals[$i][5]/($$vals[$i][5] + $$vals[$i][6] + $$vals[$i][7]) : 0);
	}
	tprint $fh, "]";
    if ( $has_vals<2 ) { return; }

	tprint $fh, "
			fig = plt.figure(figsize=($$opts{img_width},$$opts{img_height}))
			ax1 = fig.add_subplot(111)
			ax1.plot([row[0] for row in dat], [row[1] for row in dat], '-o',color='k',mec='k',markersize=3)
			ax1.set_ylabel('Number of sites',color='k')
			ax1.set_yscale('log')
			#ax1.ticklabel_format(style='sci', scilimits=(0,0), axis='y') 
			for tl in ax1.get_yticklabels(): tl.set_color('k')
			ax1.set_xlabel('Non-ref allele frequency')
			ax2 = ax1.twinx()
			ax2.plot([row[0] for row in dat], [row[2] for row in dat], '-o',color='$$opts{id2col}[$id]',mec='$$opts{id2col}[$id]',markersize=3,label='c/(c+i)')
			ax2.plot([row[0] for row in dat], [row[3] for row in dat], '--^',color='$$opts{id2col}[$id]',mec='$$opts{id2col}[$id]',markersize=3,label='c/(c+i+n/a)')
            ax2.legend(numpoints=1,markerscale=2,loc='best',prop={'size':10},frameon=False)
			ax2.set_ylabel('Indel Repeat Consistency',color='$$opts{id2col}[$id]')
			ax2.set_ylim(0,1)
			ax1.set_xlim(0,1)
			for tl in ax2.get_yticklabels(): tl.set_color('$$opts{id2col}[$id]')
			plt.subplots_adjust(right=0.88,left=0.15,bottom=0.11)
			plt.title('$$opts{title}{$id}')
			plt.savefig('$img.png')
			if pdf_plots: plt.savefig('$img.pdf')
			plt.close()
		";
}

sub plot_tstv_by_AF
{
    my ($opts,$id) = @_;
	my @vals = get_values($opts,$id,'AF');
    if ( !@vals ) { return; }

    my $fh   = $$opts{plt_fh};
	my $img  = "$$opts{prefix}tstv_by_af.$id";
	my $vals = rebin_values(\@vals,8,0);

    my $has_vals = 0;
	tprint $fh, "
			dat = [
		";
	for (my $i=0; $i<@$vals; $i++)
	{ 
		if ( $$vals[$i][2] + $$vals[$i][3] == 0 ) { next; }
        $has_vals++;
		tprint $fh, sprintf("\t[ %f, %d, %f ],\n", 
			$$vals[$i][0]/100., 
			$$vals[$i][2] + $$vals[$i][3],
			$$vals[$i][3] ? $$vals[$i][2]/$$vals[$i][3]: 0);
	}
	tprint $fh, "]";
    if ( $has_vals<2 ) { return; }

	tprint $fh, "
			fig = plt.figure(figsize=($$opts{img_width},$$opts{img_height}))
			ax1 = fig.add_subplot(111)
			ax1.plot([row[0] for row in dat], [row[1] for row in dat], '-o',color='k',mec='k',markersize=3)
			ax1.set_ylabel('Number of sites',color='k')
			ax1.set_yscale('log')
			#ax1.ticklabel_format(style='sci', scilimits=(0,0), axis='y') 
			for tl in ax1.get_yticklabels(): tl.set_color('k')
			ax1.set_xlabel('Non-ref allele frequency')
			ax2 = ax1.twinx()
			ax2.plot([row[0] for row in dat], [row[2] for row in dat], '-o',color='$$opts{id2col}[$id]',mec='$$opts{id2col}[$id]',markersize=3)
			ax2.set_ylabel('Ts/Tv',color='$$opts{id2col}[$id]')
			ax2.set_ylim(0,0.5+max(3,max(row[2] for row in dat)))
			ax1.set_xlim(0,1)
			for tl in ax2.get_yticklabels(): tl.set_color('$$opts{id2col}[$id]')
			plt.subplots_adjust(right=0.88,left=0.15,bottom=0.11)
			plt.title('$$opts{title}{$id}')
			plt.savefig('$img.png')
			if pdf_plots: plt.savefig('$img.pdf')
			plt.close()
		";
}

sub plot_tstv_by_QUAL
{
    my ($opts,$id) = @_;
	my @vals = get_values($opts,$id,'QUAL');
    if ( !@vals ) { return; }

    my $fh   = $$opts{plt_fh};
	my $img  = "$$opts{prefix}tstv_by_qual.$id";

    open(my $tfh,'>',"$img.dat") or error("$img.dat: $!");
    print $tfh "# [1]Quality\t[2]Number of sites\t[3]Marginal Ts/Tv\n";

    my @dat = ();
    my $ntot = 0;
    for my $val (@vals)
    {
        push @dat, [ $$val[0], $$val[2], $$val[3] ];    # qual, nts, ntv
        $ntot += $$val[2] + $$val[3];
    }
    my @sdat = sort { $$b[0] <=> $$a[0] } @dat;
    push @sdat, [-1];
    my $dn    = $ntot*0.05;
    my $qprev = $sdat[0][0];
    my $nts   = 0;
    my $ntv   = 0;
    my $nout  = 0;
    for my $rec (@sdat)
    {
        if ( $$rec[0]==-1 or $nts+$ntv > $dn )
        {
            if ( $ntv ) {  printf $tfh "$qprev\t%d\t%f\n", $nts+$ntv+$nout,$nts/$ntv; }
            if ( $$rec[0]==-1 ) { last; }
            $nout += $nts+$ntv;
            $nts   = 0;
            $ntv   = 0;
            $qprev = $$rec[0];
        }
        $nts += $$rec[1];
        $ntv += $$rec[2];
    }
    close($tfh) or error("close $img.dat");

	tprint $fh, "
            dat = []
            with open('$img.dat', 'rb') as f:
            \\treader = csv.reader(f, 'tab')
            \\tfor row in reader:
            \\t\\tif row[0][0] != '#': dat.append([float(x) for x in row])

			fig = plt.figure(figsize=($$opts{img_width},$$opts{img_height}))
			ax1 = fig.add_subplot(111)
			ax1.plot([row[1] for row in dat], [row[2] for row in dat], '^-', ms=3, mec='$$opts{id2col}[$id]', color='$$opts{id2col}[$id]')
			ax1.set_ylabel('Ts/Tv',fontsize=10)
			ax1.set_xlabel('Number of sites\\n(sorted by QUAL, descending)',fontsize=10)
            ax1.ticklabel_format(style='sci', scilimits=(-3,2), axis='x')
            ax1.set_ylim(min(2,min(row[2] for row in dat))-0.3,0.3+max(2.2,max(row[2] for row in dat)))

			plt.subplots_adjust(right=0.88,left=0.15,bottom=0.15)
			plt.title('$$opts{title}{$id}')
			plt.savefig('$img.png')
			if pdf_plots: plt.savefig('$img.pdf')
			plt.close()
		";
}

sub rebin_values
{
	my ($vals,$bin_size,$col,%args) = @_;
    my %avg  = exists($args{avg}) ? map {$_=>1} @{$args{avg}} : ();
	my $prev = $$vals[0][$col];
	my $iout = 0;
    my $nsum = 0;
	my (@dat,@out);
	for (my $i=0; $i<@$vals; $i++)
	{
		for (my $icol=0; $icol<@{$$vals[$i]}; $icol++)
		{
			if ( $icol==$col ) { next; }
			$dat[$icol] += $$vals[$i][$icol];
		}
        $nsum++;
		if ( $i+1<@$vals && $$vals[$i][$col] - $prev < $bin_size ) { next; }
		$dat[$col] = $prev;
		for (my $icol=0; $icol<@{$$vals[$i]}; $icol++) 
        { 
            $out[$iout][$icol] = $dat[$icol] ? $dat[$icol] : 0; 
            if ( $avg{$icol} && $nsum ) { $out[$iout][$icol] /= $nsum; }
        }
        $nsum = 0;
		@dat = ();
		$iout++;
		$prev = $$vals[$i][$col];
	}
	return \@out;
}

sub plot_concordance_by_AF
{
	my ($opts) = @_;
	my @vals = get_values($opts,2,'GCsAF');
	if ( !@vals ) { return; }
    
    # create a local copy and prepare r2 for rebinning
    @vals = @{ dclone(\@vals) };
    for (my $i=0; $i<@vals; $i++) { $vals[$i][7] *= $vals[$i][8]; }
	my $vals = rebin_values(\@vals,1,0);
	my $fh   = $$opts{plt_fh};
	my $img  = "$$opts{prefix}gts_by_af";
	my $img2 = "$$opts{prefix}r2_by_af";

	tprint $fh, "
			dat = [
		";
	for (my $i=0; $i<@$vals; $i++)
	{ 
		tprint $fh, sprintf("\t[ %f,%f,%f,%f,%f,%f,%f,%f,%d ],\n", 
            $$vals[$i][0]/100.,
            $$vals[$i][1]+$$vals[$i][4] ? $$vals[$i][1]/($$vals[$i][1]+$$vals[$i][4]) : 1,
            $$vals[$i][2]+$$vals[$i][5] ? $$vals[$i][2]/($$vals[$i][2]+$$vals[$i][5]) : 1,
            $$vals[$i][3]+$$vals[$i][6] ? $$vals[$i][3]/($$vals[$i][3]+$$vals[$i][6]) : 1,
            $$vals[$i][1]+$$vals[$i][4],
            $$vals[$i][2]+$$vals[$i][5],
            $$vals[$i][3]+$$vals[$i][6],
            $$vals[$i][8] ? $$vals[$i][7]/$$vals[$i][8] : 1,
            $$vals[$i][8]
            );
	}
	tprint $fh, "]

			fig = plt.figure(figsize=($$opts{img_width}*1.2,$$opts{img_height}))
			ax1 = fig.add_subplot(111)
			ax1.plot([row[0] for row in dat], [row[1] for row in dat],'.',color='$$opts{id2col}[1]',label='Hom RR')
			ax1.plot([row[0] for row in dat], [row[2] for row in dat],'.',color='$$opts{id2col}[0]',label='Het RA')
			ax1.plot([row[0] for row in dat], [row[3] for row in dat],'.',color='k',label='Hom AA')
			ax1.set_xlabel('Non-ref allele frequency')
			ax1.set_ylabel('Concordance')
			leg = ax1.legend(title='Concordance:',numpoints=1,markerscale=2,loc='best',prop={'size':9})
            leg.draw_frame(False)
            plt.setp(leg.get_title(),fontsize=9) 
			ax2 = ax1.twinx()
			ax2.plot([row[0] for row in dat], [row[4] for row in dat],color='$$opts{id2col}[1]')
			ax2.plot([row[0] for row in dat], [row[5] for row in dat],color='$$opts{id2col}[0]')
			ax2.plot([row[0] for row in dat], [row[6] for row in dat],color='k')
			ax2.set_ylabel('Number of genotypes')
			ax2.set_yscale('log')
			plt.subplots_adjust(left=0.15,right=0.83,bottom=0.11)
			plt.savefig('$img.png')
			if pdf_plots: plt.savefig('$img.pdf')
			plt.close()

            fig = plt.figure(figsize=($$opts{img_width}*1.3,$$opts{img_height}))
            ax1 = fig.add_subplot(111)
            ax2 = ax1.twinx()
            ax1.set_zorder(ax2.get_zorder()+1)
            ax1.patch.set_visible(False)
            ax2.plot([row[0] for row in dat], [row[8] for row in dat], '-o', color='r',mec='r',markersize=3)
            ax1.plot([row[0] for row in dat], [row[7] for row in dat], '-^', color='k',markersize=3)
            for tl in ax2.get_yticklabels(): tl.set_color('r')
            ax2.set_ylabel('Number of sites', color='r')
            ax2.set_yscale('log')
            ax1.set_ylabel('Mean allelic R\$^2\$', color='k')
            ax1.set_xlabel('Non-ref allele frequency')
			plt.subplots_adjust(left=0.19,right=0.83,bottom=0.11)
            plt.savefig('$img2.png')
            if pdf_plots: plt.savefig('$img2.pdf')
            plt.close()
		";
}

sub plot_concordance_by_sample
{
	my ($opts) = @_;
	my @vals = get_values($opts,2,'GCsS');
	if ( !@vals ) { return; }

	my $fh   = $$opts{plt_fh};
	my $img  = "$$opts{prefix}gts_by_sample";

	my @counts;
	tprint $fh, "
			dat = [
		";
	for (my $i=0; $i<@vals; $i++)
	{ 
		tprint $fh, sprintf("\t[ %d, %f, r'%s' ],\n", $i, $vals[$i][1], $vals[$i][0]);
	}
	tprint $fh, "]

			fig = plt.figure(figsize=(2*$$opts{img_width},$$opts{img_height}*0.7))
			ax1 = fig.add_subplot(111)
			ax1.plot([row[0] for row in dat], [row[1] for row in dat],'.',color='orange')
			ax1.set_ylabel('Non-ref discordance')
			ax1.set_ylim(0,)
            if sample_names: 
            \\t     plt.xticks([row[0] for row in dat],[row[2] for row in dat],**sample_font)
            \\t     plt.subplots_adjust(**sample_margins)
            else:
            \\t     plt.subplots_adjust(right=0.98,left=0.07,bottom=0.17)
            \\t     ax1.set_xlabel('Sample ID')
			plt.savefig('$img.png')
			if pdf_plots: plt.savefig('$img.pdf')
			plt.close()
		";
}

sub plot_counts_by_AF
{
	my ($opts) = @_;
	plot_counts_by_AF_col($opts,1,'SNP');
	plot_counts_by_AF_col($opts,4,'Indel');
}

sub plot_counts_by_AF_col
{
	my ($opts,$col,$title) = @_;

    my $fh  = $$opts{plt_fh};
	my $img = "$$opts{prefix}counts_by_af.".lc($title)."s";
    my @has_vals;

    for my $id (file_ids($opts))
    {
        my @tmp = get_values($opts,$id,'AF');
        my $vals = rebin_values(\@tmp,1,0);
        tprint $fh, "
            dat$id = [
            ";
        for my $val (@$vals) 
        { 
            if ( !$$val[$col] ) { next; }
            tprint $fh, "[ $$val[0], $$val[$col] ],\n"; 
            $has_vals[$id]++;
        }
        tprint $fh, "]\n";
    }
	tprint $fh, "
			fig = plt.figure(figsize=(2*$$opts{img_width},$$opts{img_height}*0.7))
			ax1 = fig.add_subplot(111)
			ax1.set_ylabel('Number of sites')
			ax1.ticklabel_format(style='sci', scilimits=(0,0), axis='y') 
			ax1.set_yscale('log')
			ax1.set_xlabel('Non-reference allele frequency')
            ax1.set_xlim(-0.05,1.05)
		";
    my $nfiles = 0;
    for my $id (file_ids($opts))
    {
        if ( !$has_vals[$id] or $has_vals[$id]<2 ) { next; }
		tprint $fh, "
			ax1.plot([row[0]/100. for row in dat$id], [row[1] for row in dat$id], '-o',markersize=3, color='$$opts{id2col}[$id]',mec='$$opts{id2col}[$id]',label='$$opts{title}{$id}')
		";
        $nfiles++;
    }
    if ( !$nfiles ) { return; }
	tprint $fh, "
            ax1.legend(numpoints=1,markerscale=1,loc='best',prop={'size':10},frameon=False)
			plt.title('$title count by AF')
			plt.subplots_adjust(bottom=0.2,left=0.1,right=0.95)
			plt.savefig('$img.png')
			if pdf_plots: plt.savefig('$img.pdf')
			plt.close()
		";

}

sub plot_overlap_by_AF
{
	my ($opts) = @_;
	plot_overlap_by_AF_col($opts,1,'SNP');
	plot_overlap_by_AF_col($opts,4,'Indel');
}

sub plot_overlap_by_AF_col
{
	my ($opts,$col,$title) = @_;

    my @ids  = file_ids($opts);
    if ( @ids != 3 ) { return; }

    my ($ia,$ib,$iab);
    for (my $i=0; $i<@ids; $i++)
    {
        if ( @{$$opts{dat}{ID}{$ids[$i]}[0]}>1 ) { $iab = $i; next; }
        if ( !defined $ia ) { $ia = $i; next; }
        $ib = $i;
    }

    my $fh  = $$opts{plt_fh};
	my $img = "$$opts{prefix}overlap_by_af.".lc($title)."s";
    my @has_vals;

    my @vals_a  = get_values($opts,$ia,'AF');
    my @vals_b  = get_values($opts,$ib,'AF');
    my @vals_ab = get_values($opts,$iab,'AF');

    my (%afs,%af_a,%af_ab);
    for my $val (@vals_a) { $afs{$$val[0]} = $$val[$col]; $af_a{$$val[0]} = $$val[$col]; }
    for my $val (@vals_ab) { $afs{$$val[0]} = $$val[$col]; $af_ab{$$val[0]} = $$val[$col]; }

    open(my $tfh,'>',"$img.dat") or error("$img.dat: $!");
    print $tfh "# [1]Allele frequency\t[2]Fraction of sites from $$opts{title}{$ids[$ia]} also in $$opts{title}{$ids[$ib]}\t[3]Number of sites\n";
    for my $af (sort { $a<=>$b } keys %afs)
    {
        my $a  = exists($af_a{$af})  ? $af_a{$af}  : 0; 
        my $ab = exists($af_ab{$af}) ? $af_ab{$af} : 0; 
        my $yval =  ($a+$ab) ? $ab * 100. / ($a + $ab) : 0;
        print $tfh "$af\t$yval\t" .($a+$ab). "\n";
    }
    close($tfh) or error("close $img.dat");

	tprint $fh, "
        dat = []
        with open('$img.dat', 'rb') as f:
        \\treader = csv.reader(f, 'tab')
        \\tfor row in reader:
        \\t\\tif row[0][0] != '#': dat.append(row)

        fig = plt.figure(figsize=(2*$$opts{img_width},$$opts{img_height}*0.7))
        ax1 = fig.add_subplot(111)
        ax1.plot([row[0] for row in dat], [row[1] for row in dat],'-o',markersize=3, color='$$opts{id2col}[1]',mec='$$opts{id2col}[1]')
        ax1.set_ylabel('Fraction found in $$opts{title}{$ib} [%]')
        ax1.set_xscale('log')
        ax1.set_xlabel('Non-reference allele frequency in $$opts{title}{$ia} [%]')
        ax1.set_xlim(0,11)
        ax1.set_xticks([0.1,0.2,0.5,1,2,5,10])
        ax1.set_xticklabels([0.1,0.2,0.5,1,2,5,10])
        plt.title('$title overlap by AF')
        plt.subplots_adjust(bottom=0.2,left=0.1,right=0.95)
        plt.savefig('$img.png')
        if pdf_plots: plt.savefig('$img.pdf')
        plt.close()

    ";

}


sub plot_indel_distribution
{
	my ($opts,$id) = @_;

	my @vals = get_values($opts,$id,'IDD');
	if ( !@vals ) { return; }

	# Set xlim to show 99 of indels but ignore outliers
	my @tmp;
	for my $id (file_ids($opts))
	{
		my @v = get_values($opts,$id,'IDD');
		for my $v (@v) { $tmp[ abs($$v[0]) ] += $$v[1]; }
	}
	my $n; 
	for my $t (@tmp) { $n += $t ? $t : 0; }
	my ($sum,$xlim);
	for ($xlim=0; $xlim<@tmp; $xlim++)
	{
		$sum += $tmp[$xlim] ? $tmp[$xlim] : 0;
		if ( $sum/$n >= 0.99 ) { last; }
	}
	if ( $xlim<20 ) { $xlim=20; }

    my $fh  = $$opts{plt_fh};
	my $img = "$$opts{prefix}indels.$id";

	tprint $fh, "
			dat = [
		";
	for my $val (@vals) { tprint $fh, "\t[$$val[0],$$val[1]],\n"; }
	tprint $fh, "]

			fig = plt.figure(figsize=($$opts{img_width},$$opts{img_height}))
			ax1 = fig.add_subplot(111)
			ax1.bar([row[0]-0.5 for row in dat], [row[1] for row in dat], color='$$opts{id2col}[0]')# , edgecolor='$$opts{id2col}[0]')
			ax1.set_xlabel('InDel Length')
			ax1.set_ylabel('Count')
			ax1.ticklabel_format(style='sci', scilimits=(0,0), axis='y')
			ax1.set_xlim(-$xlim,$xlim)
			plt.subplots_adjust(bottom=0.17)
			plt.title('$$opts{title}{$id}')
			plt.savefig('$img.png')
			if pdf_plots: plt.savefig('$img.pdf')
			plt.close()
		";
}

sub plot_substitutions
{
	my ($opts,$id) = @_;

	my @vals = get_values($opts,$id,'ST');
    if ( !@vals ) { return; }

    my $fh  = $$opts{plt_fh};
	my $img = "$$opts{prefix}substitutions.$id";

	tprint $fh, "
			dat = [
		";
	for (my $i=0; $i<@vals; $i++) { my $val=$vals[$i]; tprint $fh, "\t[$i,'$$val[0]',$$val[1]],\n"; }
	tprint $fh, "]

			fig = plt.figure(figsize=($$opts{img_width},$$opts{img_height}))
			cm  = mpl.cm.get_cmap('autumn')
			n = 12
			col = range(n)
			for i in range(n): col[i] = cm(1.*i/n)
			ax1 = fig.add_subplot(111)
			ax1.bar([row[0] for row in dat], [row[2] for row in dat], color=col)
			ax1.set_ylabel('Count')
			ax1.ticklabel_format(style='sci', scilimits=(0,0), axis='y')
			ax1.set_xlim(-0.5,n+0.5)
			plt.xticks([row[0] for row in dat],[row[1] for row in dat],rotation=45)
			plt.title('$$opts{title}{$id}')
			plt.savefig('$img.png')
			if pdf_plots: plt.savefig('$img.pdf')
			plt.close()
		";
}

sub singletons
{
	my ($opts,$id) = @_;
	my @si_vals   = get_values($opts,$id,'SiS');
	my $si_snps   = $si_vals[0][1];
	my $si_indels = $si_vals[0][4];
	my $si_irc    = sprintf "%.3f", $si_vals[0][6] ? $si_vals[0][5]/($si_vals[0][5]+$si_vals[0][6]) : 0;
	my $si_tstv   = sprintf "%.2f", $si_vals[0][3] ? $si_vals[0][2]/$si_vals[0][3] : 0;
	my @all_vals  = get_values($opts,$id,'AF');
	my $nsnps = 0;
	my $nindels = 0;
	for my $val (@all_vals)
	{
		$nsnps += $$val[1];
		$nindels += $$val[4];
	}
	$si_snps   = sprintf "%.1f", $nsnps ? $si_snps*100./$nsnps : 0;
	$si_indels = sprintf "%.1f", $nindels ? $si_indels*100./$nindels : 0;
	return { snps=>$si_snps, indels=>$si_indels, tstv=>$si_tstv, irc=>$si_irc };
}

sub calc_3n_n3n
{
    my (@vals) = @_;
    my $n3  = 0;
    my $nn3 = 0;
    for my $val (@vals)
    {
        if ( !($$val[0]%3) ) { $n3++; }
        else { $nn3++; }
    }
    if ( !$nn3 ) { return '-'; }
    return sprintf("%.2f", $n3/$nn3);
}

sub fmt_slide3v
{
	my ($opts, $image, $title) = @_;

	my $n = 0;
	for my $id (0..2) 
	{ 
		if ( -e "$image.$id.$$opts{fmt}" ) { $n++; } 
	}
	if ( !$n ) { return ''; }
	my $h = $$opts{tex}{slide3v}{"height$n"};
	my $slide = q[\vbox{];
	for my $id (0..2) 
	{ 
		if ( !-e "$image.$id.$$opts{fmt}" ) { next; }
		$slide .= qq[\\centerline{\\includegraphics[$$opts{ext},height=$h]{$image.$id}}];
	}
	$slide .= '}';
	return qq[
			% $title
			%
			\\hslide{$title}{$slide}
		];
}
sub fmt_slide3h
{
	my ($opts, $image, $title) = @_;
	my $n = 0;
	for my $id (0..2) 
	{ 
		if ( -e "$image.$id.$$opts{fmt}" ) { $n++; } 
	}
	if ( !$n ) { return ''; }
	my $w = $$opts{tex}{slide3h}{"width$n"};
	my $slide = '';
	for my $id (0..2) 
	{ 
		if ( !-e "$image.$id.$$opts{fmt}" ) { next; }
		$slide .= qq[\\includegraphics[$$opts{ext},width=$w]{$image.$id}];
	}
	return qq[
			% $title
			%
			\\hslide{$title}{$slide}
		];
}

sub bignum
{
	my ($num) = @_;
    if ( !defined $num ) { return ''; }
	if ( !($num=~/^\d+$/) ) { return $num; }
	my $len = length($num);
	my $out;
	for (my $i=0; $i<$len; $i++)
	{
		$out .= substr($num,$i,1);
		if ( $i+1<$len && !(($len-$i-1)%3) ) { $out .= ','; }
	}
	return $out;
}

sub create_pdf
{
	my ($opts) = @_;

	chdir($$opts{dir});

	my @ids     = file_ids($opts);
	my $width   = "25.4cm";	# todo: move all this to $$opts{tex}
	my $height  = "19cm";
	my $height1 = "13cm";
	my $width1  = "23cm";
	my $width2  = @ids==3 ? "10.5cm" : "10.5cm";
	my $width3  = @ids==3 ? "8cm" : "15cm";
    my $fmt     = $$opts{rasterize} ? 'png' : 'pdf';
	my $ext     = "type=$fmt,ext=.$fmt,read=.$fmt";
	my $args    = { ext=>$ext, width3=>$width3, n=>scalar @ids };
    $$opts{fmt} = $fmt;
    $$opts{ext} = $ext;

	# Check that xcolor is available
	my @has_xcolor = `kpsewhich xcolor.sty`;
	if ( !@has_xcolor )
	{
		warn("Note: The xcolor.sty package not available, black and white tables only...\n\n");
	}

	my $tex_file = "$$opts{lprefix}summary.tex";
	my $pdf_file = "$$opts{prefix}summary.pdf";
	open(my $tex,'>',$tex_file) or error("$tex_file: $!");
	tprint $tex, qq[
			% This file was produced by plot-vcfstats, the command line was:
			%   $$opts{args}
			%
			% Edit as necessary and recreate the PDF by running
			%	pdflatex $tex_file
			%

			% Slides style and dimensions
			%
			\\nonstopmode
			\\documentclass[17pt]{memoir}
			\\setstocksize{$height}{$width}
			\\settrimmedsize{\\stockheight}{\\stockwidth}{*}
			\\settrims{0pt}{0pt}
			\\setlrmarginsandblock{1cm}{*}{*}
			\\setulmarginsandblock{1.5cm}{*}{*}
			\\setheadfoot{1mm}{1cm}
			\\setlength{\\parskip}{0pt}
			\\setheaderspaces{*}{1mm}{*}
			\\setmarginnotes{1mm}{1mm}{1mm}
			\\checkandfixthelayout[fixed]
			\\usepackage{charter}	% font
			\\pagestyle{plain}
			\\makeevenfoot{plain}{}{}{\\thepage}
			\\makeoddfoot{plain}{}{}{\\thepage}
			\\usepackage{graphicx}

			% For colored tables. If xcolor.sty is not available on your system,
			% download xcolor.sty.gz LaTeX class style from 
			%	http://www.ukern.de/tex/xcolor.html
			% Unpack and install system-wide or place elsewhere and make available by
			% setting the TEXINPUTS environment variable (note the colon)
			%	export TEXINPUTS=some/dir:
			% The list of the recognised path can be obtained by running `kpsepath tex`
			%
			\\usepackage{multirow}
			\\setlength{\\tabcolsep}{0.6em}
			\\renewcommand{\\arraystretch}{1.2}
		];
	if ( @has_xcolor )
	{
		tprint $tex, '\usepackage[table]{xcolor}';
	}
	else
	{
		tprint $tex, qq[
			\\newcommand{\\definecolor}[3]{}
			\\newcommand{\\columncolor}[1]{}
			\\newcommand{\\rowcolors}[4]{}
			\\newcommand{\\arrayrulecolor}[1]{}
		];
	}
	tprint $tex, qq[
			\\definecolor{hcol1}{rgb}{1,0.6,0}
			\\definecolor{hcol2}{rgb}{1,0.68,0.2}
			\\definecolor{row1}{rgb}{1,0.88,0.7}
			\\definecolor{row2}{rgb}{1,0.92,0.8}	% #FFEBCC
			\\setlength{\\arrayrulewidth}{1.5pt}

			% Slide headings
			\\newcommand*{\\head}[1]{{\\Large\\centerline{#1}\\vskip0.5em}}

			% Slide definition
			\\newcommand*{\\hslide}[2]{%
					\\head{#1}%
					\\begin{vplace}[0.5]\\centerline{#2}\\end{vplace}\\newpage}
			\\newcommand{\\pdf}[2]{\\IfFileExists{#2.$fmt}{\\includegraphics[#1]{#2}}{}}


			% The actual slides
			\\begin{document}
		];


	# Table with summary numbers
	my $slide .= q[
            \begin{minipage}{\textwidth}\centering
            \small \rowcolors*{3}{row2}{row1} \arrayrulecolor{black}
            \begin{tabular}{l | r r r | r r r | r | r}
            \multicolumn{1}{>{\columncolor{hcol1}}l|}{}
            & \multicolumn{3}{>{\columncolor{hcol1}}c|}{SNPs}
            & \multicolumn{3}{>{\columncolor{hcol1}}c|}{indels}
            & \multicolumn{1}{>{\columncolor{hcol1}}c|}{MNPs}
            & \multicolumn{1}{>{\columncolor{hcol1}}c}{others}  \\\\
            %
            \multicolumn{1}{>{\columncolor{hcol2}}l|}{Callset}
            & \multicolumn{1}{>{\columncolor{hcol2}}c}{n}
            & \multicolumn{1}{>{\columncolor{hcol2}}c }{ts/tv}
            & \multicolumn{1}{>{\columncolor{hcol2}}c|}{\\footnotesize(1st ALT)}
            & \multicolumn{1}{>{\columncolor{hcol2}}c}{n}
            & \multicolumn{1}{>{\columncolor{hcol2}}c}{frm$^*$}
            & \multicolumn{1}{>{\columncolor{hcol2}}c|}{rc$^{**}$}
            & \multicolumn{1}{>{\columncolor{hcol2}}c|}{}
            & \multicolumn{1}{>{\columncolor{hcol2}}c}{} \\\\ \hline 
		];
    my %tex_titles;
	for my $id (@ids)
	{
		my $snps   = get_value($opts,$id,'number of SNPs:');
		my $indels = get_value($opts,$id,'number of indels:');
		my $mnps   = get_value($opts,$id,'number of MNPs:');
		my $others = get_value($opts,$id,'number of others:');
		my $tstv   = sprintf "%.2f",get_values($opts,$id,'TSTV',0,2);
		my $tstv1  = sprintf "%.2f",get_values($opts,$id,'TSTV',0,5);
		my @frsh   = get_values($opts,$id,'FS');
		my $frsh   = @frsh ? $frsh[0][3] : '--';
        my @rc     = get_values($opts,$id,'ICS');
        my $rc     = @rc ? sprintf("%.3f",$rc[0][3]) : '--';
		my $title  = $$opts{title}{$id};
		$title =~ s/_/\\_/g;
        $title =~ s/^\s*\*/\$*\$/;    # leading asterisks is eaten by TeX
        $tex_titles{$id} = $title;
		$slide .= qq[ $title & ] . bignum($snps) . qq[ & $tstv & $tstv1 & ] . bignum($indels) . qq[ & $frsh & $rc & ] . bignum($mnps) . ' & ' . bignum($others) . qq[ \\\\ \n];
	}
    $slide .= q[%
        \multicolumn{8}{r}{$^*$ frameshift ratio: out/(out+in);\hspace{1em} $^{**}$ repeat-consistency} \\\\
        \end{tabular}
        \\\\ \vspace{1em}
        \begin{tabular}{l | r r r r | r r}
        \multicolumn{1}{>{\columncolor{hcol1}}l|}{}
        & \multicolumn{4}{>{\columncolor{hcol1}}c|}{singletons {\footnotesize(AC=1)}}
        & \multicolumn{2}{>{\columncolor{hcol1}}c}{multiallelic}  \\\\
        %
        \multicolumn{1}{>{\columncolor{hcol2}}l|}{Callset}
        & \multicolumn{1}{>{\columncolor{hcol2}}c}{SNPs}
        & \multicolumn{1}{>{\columncolor{hcol2}}c}{ts/tv} 
        & \multicolumn{1}{>{\columncolor{hcol2}}c}{indels}
        & \multicolumn{1}{>{\columncolor{hcol2}}c|}{rc}
        & \multicolumn{1}{>{\columncolor{hcol2}}c}{sites}
        & \multicolumn{1}{>{\columncolor{hcol2}}c}{SNPs} \\\\ \hline 
    ];
    for my $id (@ids)
    {
        my $snps  = get_value($opts,$id,'number of SNPs:');
        my $s     = singletons($opts,$id);
        my $mals  = get_value($opts,$id,'number of multiallelic sites:');
        my $msnps = get_value($opts,$id,'number of multiallelic SNP sites:');
		my $title = $tex_titles{$id};
		$slide .= qq[ $title & $$s{snps}\\% & $$s{tstv} & $$s{indels}\\% & $$s{irc} &  ] . bignum($mals) . ' &' . bignum($msnps) . qq[ \\\\ \n];
    }
	$slide .= q[ \\end{tabular}
            \\vspace{2em}
            \\begin{itemize}[-]
            \\setlength{\\itemsep}{0pt}
        ];
    for my $id (@ids)
    {
        my $fname = $$opts{dat}{ID}{$id}[0][0];
        if ( $$opts{title}{$id} =~ / \+ / ) { next; }
        $fname =~ s/.{80}/$&\\\\\\hskip2em /g;
        $fname =~ s/_/\\_/g;
        $slide .= qq[\\item $tex_titles{$id} .. \\texttt{\\footnotesize $fname}\n];
    }
    $slide .= q[\\end{itemize}\\end{minipage}];
    my $title = exists($$opts{main_title}) ? $$opts{main_title} : 'Summary Numbers';
	tprint $tex, qq[
			% Table with summary numbers
			%
			\\hslide{$title}{$slide}
		
		];


	# Venn bars
	if ( @ids==3 )
	{
		tprint $tex, qq[%

			% Venn numbers
			%
			\\hslide{Total counts}{%
				\\includegraphics[$ext,width=$width2]{$$opts{lprefix}venn_bars.snps}%
				\\includegraphics[$ext,width=$width2]{$$opts{lprefix}venn_bars.indels}
			}
		];
	}

	tprint $tex, fmt_slide3v($opts, "$$opts{lprefix}tstv_by_sample", 'Ts/Tv by sample');
	tprint $tex, fmt_slide3v($opts, "$$opts{lprefix}hets_by_sample", 'Hets vs non-ref Homs by sample');
	tprint $tex, fmt_slide3v($opts, "$$opts{lprefix}singletons_by_sample", 'Singletons by sample {\normalsize(hets and homs)}');
	tprint $tex, fmt_slide3v($opts, "$$opts{lprefix}dp_by_sample", 'Average depth by sample');
	tprint $tex, fmt_slide3v($opts, "$$opts{lprefix}snps_by_sample", 'Number of SNPs by sample');
	tprint $tex, fmt_slide3v($opts, "$$opts{lprefix}indels_by_sample", 'Number of indels by sample');
	if ( scalar get_values($opts,2,'GCsS') ) 
	{ 
		tprint $tex, qq[
			% Genotype discordance by sample
			%
			\\hslide{Genotype discordance by sample}{\\pdf{$ext,width=$width1}{$$opts{lprefix}gts_by_sample}}
			
			]; 
	}
	if ( scalar get_values($opts,2,'GCsAF') ) 
	{
		my @vals = get_values($opts,2,'NRDs');
		my $nrd = sprintf "%.2f", $vals[0][0];
		my $rr  = sprintf "%.2f", $vals[0][1];
		my $ra  = sprintf "%.2f", $vals[0][2];
		my $aa  = sprintf "%.2f", $vals[0][3];
		my $nsamples = get_value($opts,2,'Number of samples:');
		my $table = qq[%
			{\\small
			\\rowcolors*{1}{row2}{row1}\\arrayrulecolor{black}
			\\begin{tabular}{c | c | c | c }
			\\multicolumn{1}{>{\\columncolor{hcol1}}c|}{REF/REF} & 
			\\multicolumn{1}{>{\\columncolor{hcol1}}c|}{REF/ALT} & 
			\\multicolumn{1}{>{\\columncolor{hcol1}}c|}{ALT/ALT} & 
			\\multicolumn{1}{>{\\columncolor{hcol1}}c}{NRDs} \\\\ \\hline
			$rr\\% & $ra\\% & $aa\\% & $nrd\\% \\\\
			\\end{tabular}}];
		tprint $tex, qq[
				% Genotype discordance by AF
				%
				\\head{Genotype discordance by AF}\\begin{vplace}[0.7]\\centerline{$table}%
				\\centerline{\\pdf{$ext,height=$height1}{$$opts{lprefix}gts_by_af}}\\end{vplace}
				\\newpage

                % dosage r2 by AF
                %
                \\hslide{Allelic R\$^2\$ by AF}{\\pdf{$ext,height=$height1}{$$opts{lprefix}r2_by_af}}
			];
	}
    if ( -e "$$opts{lprefix}counts_by_af.snps.$fmt" && -e "$$opts{lprefix}counts_by_af.indels.$fmt" )
    {
	    tprint $tex, qq[
	    	% SNP and indel counts by AF
	    	%
	    	\\hslide{}{\\vbox{\\noindent\\includegraphics[$ext,width=$width1]{$$opts{lprefix}counts_by_af.snps}\\\\%
	    		\\noindent\\includegraphics[$ext,width=$width1]{$$opts{lprefix}counts_by_af.indels}}}
	    	];
    }
    if ( -e "$$opts{lprefix}overlap_by_af.snps.$fmt" && -e "$$opts{lprefix}overlap_by_af.indels.$fmt" )
    {
	    tprint $tex, qq[
	    	% SNP and indel overlap by AF
	    	%
	    	\\hslide{}{\\vbox{\\noindent\\includegraphics[$ext,width=$width1]{$$opts{lprefix}overlap_by_af.snps}\\\\%
	    		\\noindent\\includegraphics[$ext,width=$width1]{$$opts{lprefix}overlap_by_af.indels}}}
	    	];
    }
	tprint $tex, fmt_slide3h($opts, "$$opts{lprefix}tstv_by_af", 'Ts/Tv by AF');
	tprint $tex, fmt_slide3h($opts, "$$opts{lprefix}tstv_by_qual", 'Ts/Tv stratified by QUAL');
	tprint $tex, fmt_slide3h($opts, "$$opts{lprefix}indels", 'Indel distribution');
	tprint $tex, fmt_slide3h($opts, "$$opts{lprefix}depth", 'Depth distribution');
	tprint $tex, fmt_slide3h($opts, "$$opts{lprefix}hwe", 'Number of HETs by AF');
	tprint $tex, fmt_slide3h($opts, "$$opts{lprefix}substitutions", 'Substitution types');
	tprint $tex, fmt_slide3h($opts, "$$opts{lprefix}irc_by_af", 'Indel Repeat Consistency by AF');
	tprint $tex, fmt_slide3h($opts, "$$opts{lprefix}irc_by_rlen", 'Indel Consistency by Repeat Type');

	tprint $tex, "\n\n\\end{document}\n";
	close($tex);

	$tex_file =~ s{^.+/}{};
	my $cmd = "pdflatex $tex_file >$$opts{logfile} 2>&1";
	print STDERR "Creating PDF: $cmd\n" unless !$$opts{verbose};
	system($cmd);
	if ( $? ) { error("The command exited with non-zero status, please consult the output of pdflatex: $$opts{dir}$$opts{logfile}\n\n"); }
	print STDERR "Finished: $pdf_file\n" unless !$$opts{verbose};
}


sub merge_vcfstats
{
	my ($opts) = @_;

    my $fh = *STDOUT;
    if ( !$$opts{merge} ) 
    { 
        open($fh,'>',"$$opts{prefix}merge.chk") or error("$$opts{prefix}merge.chk: $!\n");
    }

    print $fh "# This file was produced by plot-vcfstats, the command line was:\n#   $$opts{args}\n#\n";

    for my $sec (@{$$opts{sections}})
    {
        my $sid = $$sec{id};
        if ( !exists($$opts{dat}{$sid}) ) { next; }

        print $fh "# $$sec{header}\n$$sec{exp}\n";
        for my $id (sort keys %{$$opts{dat}{$sid}})
        {
            for my $rec (@{$$opts{dat}{$sid}{$id}})
            {
                print $fh "$sid\t$id\t", join("\t",@$rec), "\n";
            }
        }

        if ( $sid eq 'ID' )
        {
            print $fh "# $$opts{id2sec}{SN}{header}\n$$opts{id2sec}{SN}{exp}\n";
            # output summary numbers here
            for my $id (keys %{$$opts{dat}})
            {
                if ( exists($$opts{exp}{$id}) ) { next; }
                for my $key (keys %{$$opts{dat}{$id}})
                {
                    print $fh "SN\t$id\t$key\t$$opts{dat}{$id}{$key}\n";
                }
            }
        }
    }
}