#!/usr/bin/env perl
# PROJECT: CASAVA
# MODULE:  $RCSfile: generateReports.pl,v $
# AUTHOR:  Richard Carter, Lukasz Szajkowski
#
# Copyright (c) 2008 Illumina
# This software is covered by the "Illumina Genome Analyzer Software
# License Agreement" and the "Illumina Source Code License Agreement",
# and certain third party copyright/licenses, and any user of this
# source file is bound by the terms therein (see accompanying files
# Illumina_Genome_Analyzer_Software_License_Agreement.pdf and
# Illumina_Source_Code_License_Agreement.pdf and third party
# copyright/license notices).
use warnings FATAL => 'all';
use strict;
use POSIX qw(strftime);

use IO::File;
use Sys::Hostname;
use Getopt::Long;
use Carp;

use lib '/home/psgendb/local/pkg/CASAVA_v1.8.2-build/lib/CASAVA-1.8.2/perl';

use perlchartdir;

use Casava::Common::Log;
use Casava::Common::Template;
use Casava::PostAlignment::Sequencing::Config qw(loadConfiguration
  %CONF_APP %CONF_PROJ %runsConfig readProjectParameters);

use Casava::Common::IOLib
    qw(executeCmd createDirs
  file2table file2tableEx formatValueWithDeviation);
use Casava::PostAlignment::Sequencing::Reports qw(getColIndex htmlFileMenuRegister makeChart updateTotalRow);



sub buildHomHetStatsMenu($$$) {
    my ($readSampleRateStep, $readSampleRateEnd, $readSampleRateStart) = @_;

    my $ret = '';
    for (my $sampleRate = $readSampleRateStart ; $sampleRate <= $readSampleRateEnd ;
            $sampleRate += $readSampleRateStep)
    {
        my $sampleRateSetSuffix         = '';
        my $sampleRateSetName         = '';
        if ( 1.0 != $readSampleRateStart ) {
            $sampleRateSetName      = 100 * $sampleRate . 'pct';
            $sampleRateSetSuffix      = "_$sampleRateSetName";
        }

        $ret .= <<"SNIPPET_END";
<li>
    <a  id="viewPageLink"  href="snps$sampleRateSetSuffix.html"  
    class="current" accessKey="C">snps $sampleRateSetName</a>
</li>
SNIPPET_END
    }
    return $ret;
}



sub createLaneRow {
    croak "ERROR: createLaneRow " unless ( @_ == 4 );
    my ( $read, $lane, $runSummaryImg, $summaryXmlMapRef ) = @_;

    my $laneFromSummary =
      $summaryXmlMapRef->{LaneResultsSummary}{Read}[ $read - 1 ]{Lane} [ $lane - 1 ];
    my $laneParamsFromSummary =
      $summaryXmlMapRef->{LaneParameterSummary}{Lane}[ $lane - 1 ];

    my $laneReadLength = ( "unknown" eq $laneParamsFromSummary->{lengthsList} )? '' : 
    ( split( /,/, $laneParamsFromSummary->{lengthsList} ) )[ $read - 1 ];

    if (defined $laneFromSummary->{laneNumber} )
    {
      my @row = ('<p align=right>Lane ' . $laneFromSummary->{laneNumber} . " Read ${read}</p>",
                formatValueWithDeviation($laneFromSummary->{clusterCountPF} ),
                formatValueWithDeviation($laneFromSummary->{percentClustersPF} ),
                formatValueWithDeviation($laneFromSummary->{percentUniquelyAlignedPF}),
                formatValueWithDeviation( $laneFromSummary->{errorPF} ),
                $laneFromSummary->{laneYield},
                $laneReadLength,
                $runSummaryImg );
        return \@row;
    }
    else
    {
        my @row='';
        return \@row;
    }
}

sub createLaneRows($;\@;) {
    croak "ERROR: createLaneRows " unless ( @_ == 2 );
    my ( $summaryXMl, $lanes ) = @_;

    my @rows;

    my $xs = new XML::Simple( forcearray => [qw(Read Lane)] );
    my $summaryXmlMapRef = $xs->XMLin($summaryXMl);

    my $runSummaryImg =
        '<img class="emoticon" src="images/check.gif" alt="" '
      . 'align="absmiddle" border="0" height="16" width="16">';

    foreach my $lane ( @{$lanes} ) {

        my $row = createLaneRow( 1, $lane, $runSummaryImg, $summaryXmlMapRef );
        push @rows, $row;

        if ( defined $summaryXmlMapRef->{LaneResultsSummary}{Read}[1] ) {
            my $row = createLaneRow( 2, $lane, $runSummaryImg, $summaryXmlMapRef );
            push @rows, $row;
        }

    }

    return @rows;
}



my $usage =
  "generateReports.pl [options]\n" .
  "\t--projectDir - Project directory\n" .
  "\t--target     - target which part of the script to run\n"
;
my $projectDir = "";
my $help       = 0;
my $target     = 'all';
my $sampleRateSetSuffix = '';
my $result = GetOptions(
                         'projectDir|p=s' => \$projectDir,
                         'target=s'       => \$target,
                         'sampleRateSetSuffix=s' => \$sampleRateSetSuffix,
                         'help|h'         => \$help
);
errorExit "ERROR: Incorrect number of arguments\n$usage" if ( $projectDir eq "" );
loadConfiguration( $projectDir );
#############################################################
# Configuration
my $cmd               = "";
my $timeStampFormat   = $CONF_APP{formatTimeStamp};
my $hostname          = hostname();
my $argvStr           = join ' ', @ARGV;
my $projectResultsDir = $CONF_PROJ{dirBuildExport};
my $buildResultsDir   = $CONF_PROJ{dirBuildParsed};
my $readSampleRateStart    = $CONF_PROJ{readSampleRateStart};
my $readSampleRateEnd      = $CONF_PROJ{readSampleRateEnd};
my $readSampleRateStep     = $CONF_PROJ{readSampleRateStep};

my $readMode        = $CONF_PROJ{readMode};
my $readCode        = ( $readMode eq 'single' ? "SE" : "PE");
my $projectName = "CASAVA $CONF_PROJ{appVersion} $readCode " . $CONF_PROJ{applicationType} . "-Seq ";

#############################################################
if ($help) {
    print $usage;
    exit(0);
}
printLog(
          "Running $hostname:[$0 " . $argvStr . "] at "
            . ( strftime $timeStampFormat, localtime ) . "\n",
          5
);
my $templPath = File::Spec->catdir( '/home/psgendb/local/pkg/CASAVA_v1.8.2-build/share/CASAVA-1.8.2', "templates" );
my $htmlPath  = File::Spec->catdir( $projectDir,         "html" );
my $statsPath = File::Spec->catdir( $projectDir,         "stats" );
configure($templPath);
my %projectChromsBinSizes = ();
readProjectParameters( %projectChromsBinSizes, "PROJECT_BIN_SIZES", $projectDir );
my @chromsProject = sort keys %projectChromsBinSizes;

if ( $target eq 'configure' || $target eq 'all' ) {
    my @dirs = ( 'html', 'stats' );
    createDirs( "$projectDir", @dirs );
    my @cpdirs = qw(css js images);
    # simple copy from the installation location sometimes creates files that are group-read-only
    # evne if umask is 002. Add all-write during tarring and make untarring apply umask to make the
    # copied files look like the ones that get created by the build (fix for BF-901)
    my $cmd = "tar --mode='a+w' -cO -C $templPath @cpdirs |tar --no-same-permissions -x -C $htmlPath";
    executeCmd( $cmd, 4 );
    my $title       = "$projectName - analysis";
    my %values      = ();
    my $path        =
        "<a href=\"Home.html\">$projectName</a>"
      . "&gt;&nbsp;<a href=\"copyright.html\">copyright</a>";
    $values{path}        = $path;
    $values{title}       = $title;
    $values{projectName} = $projectName;
    my $filename = File::Spec->catfile( $htmlPath, 'copyright.html' );
    my $html = fromHtmlTemplate2HtmlFile( %values, 'copyright', $filename );
}



if ( $target eq 'main' || $target eq 'all' ) {

    my $title       = "$projectName - analysis";
    my %values      = ();
    my $path        = "<a href=\"Home.html\">$projectName</a>";
    $values{path}        = $path;
    $values{title}       = $title;
    $values{projectName} = $projectName;
    my @table = ();

    # check all the files exist
    for my $runIdTmp ( sort keys %{ $runsConfig{run} } ) {
        my %sets = %{ $runsConfig{run}->{$runIdTmp}->{set} };
        foreach my $setIdTmp ( keys %sets ) {
            my @lanes =
              @{ $runsConfig{run}->{$runIdTmp}->{set}->{$setIdTmp}->{lanes} };
            my $geraldFullDir =
              $runsConfig{run}->{$runIdTmp}->{set}->{$setIdTmp}->{gerald};
            my $lanesStr = join( ',&nbsp;', @lanes );
            my $runSummaryFile =
              File::Spec->catfile( $geraldFullDir, 'Summary.htm' );
            my $runSummary    = '';
            my $runSummaryImg = '';
            if ( -e $runSummaryFile ) {
                $runSummary =
                    "<a alt=\"$geraldFullDir\" "
                  . "href=\"$geraldFullDir/Summary.htm\">$runIdTmp</a>";
                $runSummary =~ s/\/lustre\//http\:\/\//g;
                $runSummaryImg =
                    '<img class="emoticon" src="images/check.gif" alt="" '
                  . 'align="absmiddle" border="0" height="16" width="16">';
            }
            else {
                $runSummary = $runIdTmp;
            }

            my $runSummaryXmlFile =
              File::Spec->catfile( $geraldFullDir, 'Summary.xml' );

            if ( -e $runSummaryXmlFile ) {
                if ( 0 == @table ) {
                    my @header = (
                                   'Name',
                                   'Clusters (PF)',
                                   '%PF Clusters',
                                   '%Align (PF)',
                                   '%Error Rate (PF)',
                                   'Yield',
                                   'Read Length',
                                   'Status'
                    );
                    push @table, \@header;
                }

                my @row = ( $runSummary, '', '', '', '', '', '', '' );
                push @table, \@row;

                my @rows = createLaneRows( $runSummaryXmlFile, @lanes );
                push @table, @rows;
            }
            else {
                if ( 0 == @table ) {
                    my @header = (
                                   'Run name',
                                   'cycles  columns',
                                   'PF yield (whole FC)',
                                   'Used Yield (Gb)',
                                   'use lanes',
                                   'Summary'
                    );
                    push @table, \@header;
                }

                my @row = 
                    ( $runSummary, 'N/A', 'N/A', 'N/A', $lanesStr, $runSummaryImg );
                push @table, \@row;
            }
        }
    }    # for
    my $tableHtml = htmlTable( @table, 'confluence' );
    $values{table} = $tableHtml;

    my $htmlFileName = 'Home.html';
    my $htmlFilePath = File::Spec->catfile( $htmlPath, $htmlFileName );

    $values{page} = $htmlFileName;
    my $html = fromHtmlTemplate2HtmlFile( %values, 'body', $htmlFilePath );
    htmlFileMenuRegister( $htmlPath , $htmlFileName, "Sequence");
}



sub sprintRatio($$) {
    my ($n,$d) = @_;
    return sprintf( "%5.5f", ($d==0 ? 0. : ($n/$d)));
}



sub formatDuplicatesRow {
    my ($rowRef,$colIndexRef) = @_;

    my $nondup = $$rowRef[ $$colIndexRef{nonduplicate_pairs} ];
    my $dup = $$rowRef[ $$colIndexRef{duplicate_pairs} ];
    my $total = $nondup+$dup;
    my $frac_nondup = sprintRatio($nondup,$total);
    my $frac_dup = sprintRatio($dup,$total);

    return ($$rowRef[ $$colIndexRef{seq_name} ],
            $total,
            $nondup,
            $frac_nondup,
            $dup,
            $frac_dup);
}



if ( $target eq 'duplicates' || $target eq 'all' ) {
    my @tableHeader = (
                       'name',
                       'total read pairs',
                       'good read pairs',
                       'fraction good read pairs',
                       'duplicate read pairs',
                       'fraction duplicate read pairs'
                      );

    my @tableHtml = ( \@tableHeader );

    my $dupCountFile = File::Spec->catfile( $statsPath, 'dupCount.summary.txt' );
    my $tableHtmlStr = '';
    if ( -e $dupCountFile ) {
        my @Xaxis              = ();
        my @goodReads          = ();
        my @goodReadsPr        = ();
        my @duplicates         = ();
        my @duplicatesPr       = ();
        my @totalRow           = ();

        my @table = ();
        my %fileMD = ();
        file2tableEx(@table, %fileMD,  $dupCountFile);

        my %colIndex = ();
        getColIndex(%colIndex, %fileMD, $dupCountFile);

        for my $rowRef (@table) {
            updateTotalRow(\@totalRow,$rowRef,$colIndex{seq_name});
            my @tableRow = formatDuplicatesRow($rowRef,\%colIndex);
            push @tableHtml, \@tableRow;

            push @Xaxis,         $tableRow[ 0 ];
            push @goodReads,     $tableRow[ 2 ];
            push @goodReadsPr,   $tableRow[ 3 ];
            push @duplicates,    $tableRow[ 4 ];
            push @duplicatesPr,  $tableRow[ 5 ];
        }

        my @totalTableRow = formatDuplicatesRow(\@totalRow,\%colIndex);

        {
            my $plotTitle = "Number of read pairs by duplicate type";
            my $pngFile = "duplicates.png";
            my $xAxisLabel = "sequence name";
            my $yAxisLabel = "number of read pairs";

            makeChart($htmlPath,
                      $pngFile,
                      \$tableHtmlStr,
                      $plotTitle,
                      $xAxisLabel,
                      $yAxisLabel,
                      \@Xaxis,
                      [ [ \@goodReads,   0x80ff80, "Good Read Pairs" ],
                        [ \@duplicates,  0x8080ff, "Duplicates" ]]);
        }

        {
            my $plotTitle = "Fraction of read pairs by duplicate type";
            my $pngFile = "duplicatesPr.png";
            my $xAxisLabel = "sequence name";
            my $yAxisLabel = "fraction of read pairs";

            makeChart($htmlPath,
                      $pngFile,
                      \$tableHtmlStr,
                      $plotTitle,
                      $xAxisLabel,
                      $yAxisLabel,
                      \@Xaxis,
                      [ [ \@goodReadsPr,   0x80ff80, "Good Read Pairs" ],
                        [ \@duplicatesPr,  0x8080ff, "Duplicates" ]]);
        }

        $tableHtmlStr .= htmlTableWithTotals( @tableHtml, @totalTableRow, 'confluence' );
    }
    else {
        $tableHtmlStr = "<b>Not available</b>";
    }
    my $htmlFileName = 'dupCount.html';
    my $htmlFilePath = File::Spec->catfile( $htmlPath, $htmlFileName );
    my %values      = ();
    $values{path}      =
      "<a href=\"Home.html\">$projectName</a>"
        . "&gt;&nbsp;<a href=\"$htmlFileName\">duplicates</a>";
    $values{title}       = "$projectName - analysis: duplicates";
    $values{projectName} = $projectName;
    $values{table} = $tableHtmlStr;
    $values{page} = $htmlFileName;
    fromHtmlTemplate2HtmlFile( %values, 'body', $htmlFilePath );
    htmlFileMenuRegister( $htmlPath , $htmlFileName, "Duplicates");
}



if ( $target eq 'depth' || $target eq 'all' ) {
    for (my $sampleRate = $readSampleRateStart ; $sampleRate <= $readSampleRateEnd ;
            $sampleRate += $readSampleRateStep)
    {
        my $tableHtmlStr = '';

        my $sampleRateSetSuffix         = '';
        my $sampleRateSetName         = '';
        if ( 1.0 != $readSampleRateStart ) {
            $sampleRateSetName      = 100 * $sampleRate . 'pct';
            $sampleRateSetSuffix      = "_$sampleRateSetName";
        }


        #coverage.chrom.txt
        my $fileInSnp   = File::Spec->catfile( $statsPath, "snpChromStats$sampleRateSetSuffix.txt" );
        my $fileInDepth = File::Spec->catfile( $statsPath, "coverage.chrom$sampleRateSetSuffix.txt" );
        my @tableInSnp  = ();
        my @tableInDepth = ();
        my @tableHtml    = ();
        my @headerChrom = (
                            'Name',
                            'homozygous',
                            'heterozygous',
                            'het1',
                            'het2',
                            'other het',
                            'total SNPs',
                            'het %',
                            'total used bases',
                            'mean depth',
                            'corrected mean depth',
            'covered %',
                            'length',
                            'corrected length'
        );
        my %snpDepFields = (
            setName           => 0,
            het1Count         => 1,
            het2Count         => 2,
            otherHetCount     => 3,
            diffCount         => 4,
            setName2          => 5,
            positionsCovered  => 6,
            coverage          => 7,
            chromosomeLength  => 8,
            knownBases        => 9,
        );
        
        my %depthStatsColumns =(
                    setName             => 0,
                    homozygous          => 1,
                    heterozygous        => 2,
                    het1Count           => 3,
                    het2Count           => 4,
                    otherHetCount       => 5,
                    totalSNPs           => 6,
                    hetPr               => 7,
                    coverage            => 8,
                    depthPr             => 9,
                    knownBasesPr        => 10,
                    coveredPr           => 11,
                    chromosomeLength    => 12,
                    knownBases          => 13,
        );
        
        my $chromosomeLength = 0;
        my $knownBases       = 0;
        my $positionsCoveredTotal = 0;
        push @tableHtml, \@headerChrom;
    
        if ( -e $fileInSnp && -e $fileInDepth ) {
            my @Xaxis                  = ();
            my @Yaxis                  = ();
            my @YaxisDepth             = ();
            my @heterozygousData       = ();
            my @homozygousData         = ();
            my @meanDepthData          = ();
            my @correctedMeanDepthData = ();
            my @totalsRow = ( 'Totals', 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 );
            file2table( @tableInSnp,   $fileInSnp );
            file2table( @tableInDepth, $fileInDepth );
            my $rowId = 0;
    
            foreach my $rowRef (@tableInSnp) {
                my @row  = @{$rowRef};
                my @row2 = @{ shift @tableInDepth };
                push @row, @row2;
                my $str = join( "\t", @row );
                my $homozygous = $row[ $snpDepFields{diffCount} ];
                my $heterozygous = $row[ $snpDepFields{het1Count} ] +
                  $row[ $snpDepFields{het2Count} ] +
                  $row[ $snpDepFields{otherHetCount} ];
                my $totalSNP = $homozygous + $heterozygous;
                my $hetPr;
    
                $hetPr = sprintf( "%5.2f%%", $totalSNP ? ($heterozygous / $totalSNP * 100) : 0 );
                my $depthX = sprintf( "%5.1f",
                                      $row[ $snpDepFields{coverage} ] /
                                        $row[ $snpDepFields{chromosomeLength} ] );
    
                $chromosomeLength = $row[ $snpDepFields{chromosomeLength} ];
                $knownBases       = $row[ $snpDepFields{knownBases} ];
                my $knownBasesPr = sprintf( "%5.1f", $knownBases ?
                    ($row[ $snpDepFields{coverage} ] / $knownBases) : 0 );
    
                my $positionsCovered = $row[ $snpDepFields{positionsCovered} ];
                my $positionsCoveredPr = sprintf( "%5.1f%%", $knownBases ? 
                    ($positionsCovered / $knownBases * 100) : 0.0 );
    
                push @Xaxis,                  $row[ $snpDepFields{setName} ];
                push @Yaxis,                  $homozygous;
                push @heterozygousData,       $heterozygous;
                push @homozygousData,         $homozygous;
                push @meanDepthData,          $depthX;
                push @correctedMeanDepthData, $knownBasesPr;
    
                my @newRow = (
                    $row[ $snpDepFields{setName} ],
                    $homozygous,
                    $heterozygous,
                    $row[ $snpDepFields{het1Count} ],
                    $row[ $snpDepFields{het2Count} ],
                    $row[ $snpDepFields{otherHetCount} ],
                    $totalSNP,
                    $hetPr,
                    $row[ $snpDepFields{coverage} ],
                    $depthX,
                    $knownBasesPr,
                    $positionsCoveredPr,
                    $chromosomeLength,
                    $row[ $snpDepFields{knownBases} ]
                );
                push @tableHtml, \@newRow;
    
                $positionsCoveredTotal += $positionsCovered;
                foreach my $k (keys(%depthStatsColumns)) {
                    $totalsRow[$depthStatsColumns{$k}] += $newRow[$depthStatsColumns{$k}]
                       unless (grep (/^$k$/, ('setName', 'hetPr', 'depthPr', 'coveredPr', 'knownBasesPr')));
                }
                $rowId++;
            }
    
            $totalsRow[$depthStatsColumns{hetPr}] =  sprintf( "%5.2f%%", 
                ( $totalsRow[$depthStatsColumns{totalSNPs}] ) 
                    ? $totalsRow[$depthStatsColumns{heterozygous}] / $totalsRow[$depthStatsColumns{totalSNPs}] * 100
                    : 0 );
    
            $totalsRow[$depthStatsColumns{depthPr}] = sprintf( "%5.1f", $totalsRow[$depthStatsColumns{chromosomeLength}] ? $totalsRow[$depthStatsColumns{coverage}] / $totalsRow[$depthStatsColumns{chromosomeLength}] : 0.0 );
            my $totalCorrectedMeanDepth = $totalsRow[$depthStatsColumns{knownBases}] ? $totalsRow[$depthStatsColumns{coverage}] / $totalsRow[$depthStatsColumns{knownBases}] : 0.0;
            $totalsRow[$depthStatsColumns{knownBasesPr}] = sprintf( "%5.1f", $totalCorrectedMeanDepth );
            $totalsRow[$depthStatsColumns{coveredPr}] = sprintf( "%5.2f%%", $totalsRow[$depthStatsColumns{knownBases}] ? ($positionsCoveredTotal / $totalsRow[$depthStatsColumns{knownBases}] * 100) : 0 );
    
            #chrom M value is typically so high that all others will stop being
            #visible if we plot it.
            for ( my $depthChartBar = scalar(@correctedMeanDepthData) - 1 ;
                  $depthChartBar >= 0; --$depthChartBar )
            {
                # exclude from chart those that exceed the mean by too much.
                if ( $correctedMeanDepthData[$depthChartBar] > $totalCorrectedMeanDepth * 4 )
                {
                    splice( @Xaxis,                  $depthChartBar, 1 );
                    splice( @Yaxis,                  $depthChartBar, 1 );
                    splice( @heterozygousData,       $depthChartBar, 1 );
                    splice( @homozygousData,         $depthChartBar, 1 );
                    splice( @meanDepthData,          $depthChartBar, 1 );
                    splice( @correctedMeanDepthData, $depthChartBar, 1 );
                }
            }
    
            my $pngFile = "snpHomHet$sampleRateSetSuffix.png";
            $tableHtmlStr .= "<img src=\"$pngFile\" >";
            $tableHtmlStr .=
              "<a id=\"depth\"></a><h2>SNPs " . " for all chromosomes<h2>";
    
            # plot the chart
            my $pngFilePath = File::Spec->catfile( $htmlPath, $pngFile );
            my $c = new XYChart( 1024, 500 );
            $c->setPlotArea( 60, 80, 900, 400, 0xeeeeee, 0xffffff );
            my $layer = $c->addBarLayer2( $perlchartdir::Stack, 8 );
    
            # Add the data sets to the bar layer
            $c->xAxis()->setLabels( \@Xaxis );
            $layer->addDataSet( \@homozygousData,   0x80ff80, "homozygous" );
            $layer->addDataSet( \@heterozygousData, 0xff8080, "heterozygous" );
            $layer->setBorderColor( $perlchartdir::Transparent,
                                  perlchartdir::softLighting($perlchartdir::Left) );
            $c->yAxis()->setTitle( "number of SNPs", "arialbd.ttf", 10 );
            $c->xAxis()->setTitle( "Chromosomes",    "arialbd.ttf", 10 );
            $c->xAxis()->setLabelStyle( "arialbd.ttf", 8 );
            $c->yAxis()->setLabelStyle( "arialbd.ttf", 8 );
            $c->addTitle( "Number of SNPs for all chromosomes", "timesbi.ttf", 18 );
            $c->addLegend( 55, 32, 0 );
            $c->makeChart($pngFilePath);
            $pngFile = "depth$sampleRateSetSuffix.png";
            $tableHtmlStr .= "<img src=\"$pngFile\" >";
    
            # plot the chart
            $pngFilePath = File::Spec->catfile( $htmlPath, $pngFile );
            $c = new XYChart( 1024, 500 );
            $c->setPlotArea( 60, 80, 900, 400, 0xeeeeee, 0xffffff );
            $layer = $c->addBarLayer2( $perlchartdir::Stack, 8 );
    
            # Add the data sets to the bar layer
            $c->xAxis()->setLabels( \@Xaxis );
            $c->yAxis()->setLabels( \@YaxisDepth );
            $layer->addDataSet( \@meanDepthData, 0x80ff80, "Mean depth" );
    
    #$layer->addDataSet( \@correctedMeanDepthData, 0xff8080, "Corrected mean depth" );
            $layer->setBorderColor( $perlchartdir::Transparent,
                                  perlchartdir::softLighting($perlchartdir::Left) );
            $c->yAxis()->setTitle( "Depth",       "arialbd.ttf", 10 );
            $c->xAxis()->setTitle( "Chromosomes", "arialbd.ttf", 10 );
            $c->xAxis()->setLabelStyle( "arialbd.ttf", 8 );
            $c->yAxis()->setLabelStyle( "arialbd.ttf", 8 );
            $c->addTitle( "Depth for all chromosomes", "timesbi.ttf", 18 );
            $c->addLegend( 55, 32, 0 );
            $c->makeChart($pngFilePath);
            $tableHtmlStr .=
              "<a id=\"depth\"></a><h2>Depth" . " for all chromosomes<h2>";
    
            $tableHtmlStr .= htmlTableWithTotals( @tableHtml, @totalsRow, 'confluence' );
    
        }
        else {
            $tableHtmlStr = "<b>Not available</b>";
        }
        my %values      = ();
        $values{projectName} = $projectName;
    
        # check all the files exist
        my $htmlFileName =  "snpHomHet$sampleRateSetSuffix.html";
        my $htmlFilePath = File::Spec->catfile( $htmlPath, $htmlFileName );
        $values{title} = "$projectName - analysis: SNPs for Chromosomes $sampleRateSetName";
        $values{path}  =
            "<a href=\"Home.html\">$projectName</a>"
          . "&gt;&nbsp;<a href=\"snpHomHet$sampleRateSetSuffix.html\">SNPs and depth for Chromosomes</a>";
        $values{table} = $tableHtmlStr;
        $values{page} = $htmlFileName;
        my $html = fromHtmlTemplate2HtmlFile( %values, 'body', $htmlFilePath );


# Transfer this logic to use the htmlFileMenuRegister instead:
#
#        buildHomHetStatsMenu ($readSampleRateStep, 
#                              $readSampleRateEnd, 
#                              $readSampleRateStart);
#

        htmlFileMenuRegister( $htmlPath , $htmlFileName, "Ca1.7: Depth and SNPs");
    }
}
if ( $target eq 'coverageSimulation' || $target eq 'all' ) {
    my $tableHtmlStr = '';
    foreach my $chrom ( @chromsProject ) {
        my $fileIn    = File::Spec->catfile( $statsPath, "$chrom.covsim.txt" );
        my @tableIn   = ();
        my @tableHtml = ();

        #my $coverageFile = File::Spec->catfile( $statsPath, 'coverage.txt' );
        my @headerChrom = (
                            'Name',
                            'homozygous',
                            'heterozygous',
                            'het1',
                            'het2',
                            'other het',
                            'total SNPs',
                            'het %',
                            'total used bases',
                            'mean depth',
                            'corrected mean depth'
        );
        my %covSimFields = (
                             'setName'           => 0,
                             'het1Count'         => 1,
                             'het2Count'         => 2,
                             'otherHetCount'     => 3,
                             'diffCount'         => 4,
                             'numberOfPositions' => 5,
                             'coverage'          => 6,
                             'chromosomeLength'  => 7,
                             'knownBases'        => 8,
        );
        my $chromosomeLength = 0;
        my $knownBases       = 0;
        push @tableHtml, \@headerChrom;
        if ( -e $fileIn ) {
            my @Xaxis            = ();
            my @Yaxis            = ();
            my @heterozygousData = ();
            my @homozygousData   = ();
            file2table( @tableIn, $fileIn );
            my $rowId = 0;
            foreach my $rowRef (@tableIn) {
                my @row    = @{$rowRef};
                my @newRow = ();
                push @newRow, $row[ $covSimFields{setName} ];
                my $homozygous = $row[ $covSimFields{diffCount} ];
                push @newRow,         $homozygous;
                push @Yaxis,          $homozygous;
                push @homozygousData, $homozygous;
                my $heterozygous = $row[ $covSimFields{het1Count} ] +
                  $row[ $covSimFields{het2Count} ] +
                  $row[ $covSimFields{otherHetCount} ];
                push @heterozygousData, $heterozygous;
                push @newRow,           $heterozygous;
                push @newRow,           $row[ $covSimFields{het1Count} ];
                push @newRow,           $row[ $covSimFields{het2Count} ];
                push @newRow,           $row[ $covSimFields{otherHetCount} ];
                my $totalSNP = $homozygous + $heterozygous;
                push @newRow, $totalSNP;
                my $hetPr =
                  sprintf( "%5.2f%%", $heterozygous / $totalSNP * 100 );
                push @newRow, $hetPr;
                push @newRow, $row[ $covSimFields{coverage} ];
                my $coverageX = sprintf( "%5.1f",
                                    $row[ $covSimFields{coverage} ] /
                                      $row[ $covSimFields{chromosomeLength} ] );
                push @Xaxis,  $coverageX;
                push @newRow, $coverageX;
                $chromosomeLength = $row[ $covSimFields{chromosomeLength} ];
                $knownBases       = $row[ $covSimFields{knownBases} ];
                my $knownBasesPr = sprintf( "%5.1f",
                           $row[ $covSimFields{coverage} ] /
                             $row[ $covSimFields{knownBases} ] );
                push @newRow,    $knownBasesPr;
                push @tableHtml, \@newRow;
                $rowId++;
            }    # foreach
            my $pngFile = "$chrom.covsim.png";
            $tableHtmlStr .=
                "<a id=\"$chrom\"></a><h2>Depth for $chrom"
              . " (length = $chromosomeLength corrected length = $knownBases)<h2>";
            $tableHtmlStr .= htmlTable( @tableHtml, 'confluence' );
            $tableHtmlStr .= "<img src=\"$pngFile\" >";

            # plot the chart
            my $pngFilePath = File::Spec->catfile( $htmlPath, $pngFile );
            my $c = new XYChart( 1024, 600 );
            $c->setPlotArea( 60, 80, 900, 400, 0xeeeeee, 0xffffff );
            my $layer = $c->addBarLayer2( $perlchartdir::Stack, 8 );

            # Add the data sets to the bar layer
            $c->xAxis()->setLabels( \@Xaxis );
            $layer->addDataSet( \@homozygousData,   0x80ff80, "homozygous" );
            $layer->addDataSet( \@heterozygousData, 0xff8080, "heterozygous" );
            $layer->setBorderColor( $perlchartdir::Transparent,
                perlchartdir::softLighting($perlchartdir::Left));
            $c->yAxis()->setTitle( "number of SNPs", "arialbd.ttf", 10 );
            $c->xAxis()->setTitle( "mean depth",     "arialbd.ttf", 10 );
            $c->xAxis()->setLabelStyle( "arialbd.ttf", 8 );
            $c->yAxis()->setLabelStyle( "arialbd.ttf", 8 );
            $c->addTitle( "Number of SNPs vs mean depth for $chrom",
                          "timesbi.ttf", 18 );
            $c->addLegend( 55, 32, 0 );
            $c->makeChart($pngFilePath);
        }
        else {
            $tableHtmlStr = "<b>Not available</b>";
        }
    }
    my %values      = ();
    $values{projectName} = $projectName;

    # check all the files exist
    my $htmlFileName = 'coverageSimulation.html';
    my $htmlFilePath = File::Spec->catfile( $htmlPath, $htmlFileName );
    $values{title} = "$projectName - analysis: Depth for Chromosomes";
    $values{path}  = 
        "<a href=\"Home.html\">$projectName</a>"
      . "&gt;&nbsp;<a href=\"coverageSimulation.html\">Depth simulation</a>";
    $values{table} = $tableHtmlStr;
    $values{page} = $htmlFileName;
    my $html = fromHtmlTemplate2HtmlFile( %values, 'body', $htmlFilePath );
    htmlFileMenuRegister( $htmlPath , $htmlFileName, "Coverage Simulation");
}



sub formatIndelRow {
    my ($rowRef,$colIndexRef) = @_;

    my $homozygous = $$rowRef[ $$colIndexRef{hom} ];
    my $heterozygous = $$rowRef[ $$colIndexRef{het} ];
    my $insertion = $$rowRef[ $$colIndexRef{insertion} ];
    my $deletion = $$rowRef[ $$colIndexRef{deletion} ];
    my $complex = $$rowRef[ $$colIndexRef{complex} ];
    my $totalIndel = $homozygous + $heterozygous;
    my $hetHomRatio = sprintf( "%5.3f", $homozygous ? ($heterozygous / $homozygous) : 0 );

    return ($$rowRef[ $$colIndexRef{chromosome} ],
            $insertion,
            $deletion,
            $complex,
            $homozygous,
            $heterozygous,
            $hetHomRatio,
            $totalIndel);
}



if ( $target eq 'indels' || $target eq 'all' ) {
    my $tableHtmlStr = '';

    my $fileInIndel   = File::Spec->catfile( $statsPath, "indelChromStats.txt" );
    my @tableInIndel  = ();
    my @tableHtml     = ();
    my @headerChrom = (
                        'name',
                        'insertion',
                        'deletion',
                        'complex',
                        'homozygous',
                        'heterozygous',
                        'het/hom ratio',
                        'total indels'
    );
    
    push @tableHtml, \@headerChrom;

    if ( -e $fileInIndel ) {
        my @Xaxis                  = ();
        my @heterozygousData       = ();
        my @homozygousData         = ();
        my @insertionData          = ();
        my @deletionData           = ();
        my @complexData            = ();

        my %fileMD = ();
        file2tableEx( @tableInIndel, %fileMD,  $fileInIndel );

        my %colIndex = ();
        getColIndex(%colIndex, %fileMD, $fileInIndel);

        my @totalRow = ();
        foreach my $rowRef (@tableInIndel) {
            updateTotalRow(\@totalRow,$rowRef,$colIndex{chromosome});
            my @row  = @{$rowRef};

            push @Xaxis,                  $row[ $colIndex{chromosome} ];
            push @heterozygousData,       $row[ $colIndex{het} ];
            push @homozygousData,         $row[ $colIndex{hom} ];
            push @insertionData,          $row[ $colIndex{insertion} ];
            push @deletionData,           $row[ $colIndex{deletion} ];
            push @complexData,            $row[ $colIndex{complex} ];

            my @tableRow = formatIndelRow(\@row,\%colIndex);
            push @tableHtml, \@tableRow;
        }

        my @totalTableRow = formatIndelRow(\@totalRow,\%colIndex);

        {
            my $plotTitle = "Indel category by chromosome";
            my $pngFile = "indelCat.png";
            my $pngFilePath = File::Spec->catfile( $htmlPath, $pngFile );

            $tableHtmlStr .= "<img src=\"$pngFile\" >";
            $tableHtmlStr .=
                "<a id=\"indelCat\"></a><h2>$plotTitle<h2>";
            
            # plot the chart
            my $c = new XYChart( 1024, 500 );
            $c->setPlotArea( 60, 80, 900, 400, 0xeeeeee, 0xffffff );
            my $layer = $c->addBarLayer2( $perlchartdir::Stack, 8 );
            
            # Add the data sets to the bar layer
            $c->xAxis()->setLabels( \@Xaxis );
            $layer->addDataSet( \@insertionData,   0x80ff80, "insertion" );
            $layer->addDataSet( \@deletionData, 0xff8080, "deletion" );
            $layer->addDataSet( \@complexData, 0x8080ff, "complex" );
            $layer->setBorderColor( $perlchartdir::Transparent,
                                    perlchartdir::softLighting($perlchartdir::Left) );
            $c->yAxis()->setTitle( "number of indels", "arialbd.ttf", 10 );
            $c->xAxis()->setTitle( "chromosome",    "arialbd.ttf", 10 );
            $c->xAxis()->setLabelStyle( "arialbd.ttf", 8 );
            $c->yAxis()->setLabelStyle( "arialbd.ttf", 8 );
            $c->addTitle( $plotTitle, "timesbi.ttf", 18 );
            $c->addLegend( 55, 32, 0 );
            $c->makeChart($pngFilePath);
        }

        {
            my $plotTitle = "Indel call type by chromosome";
            my $pngFile = "indelHomHet.png";
            my $pngFilePath = File::Spec->catfile( $htmlPath, $pngFile );

            $tableHtmlStr .= "<img src=\"$pngFile\" >";
            $tableHtmlStr .=
                "<a id=\"indelHomHet\"></a><h2>$plotTitle<h2>";

            # plot the chart
            my $c = new XYChart( 1024, 500 );
            $c->setPlotArea( 60, 80, 900, 400, 0xeeeeee, 0xffffff );
            my $layer = $c->addBarLayer2( $perlchartdir::Stack, 8 );

            # Add the data sets to the bar layer
            $c->xAxis()->setLabels( \@Xaxis );
            $layer->addDataSet( \@homozygousData,   0x80ff80, "homozygous" );
            $layer->addDataSet( \@heterozygousData, 0xff8080, "heterozygous" );
            $layer->setBorderColor( $perlchartdir::Transparent,
                                    perlchartdir::softLighting($perlchartdir::Left) );
            $c->yAxis()->setTitle( "number of indels", "arialbd.ttf", 10 );
            $c->xAxis()->setTitle( "chromosome", "arialbd.ttf", 10 );
            $c->xAxis()->setLabelStyle( "arialbd.ttf", 8 );
            $c->yAxis()->setLabelStyle( "arialbd.ttf", 8 );
            $c->addTitle( $plotTitle, "timesbi.ttf", 18 );
            $c->addLegend( 55, 32, 0 );
            $c->makeChart($pngFilePath);
        }

        $tableHtmlStr .= htmlTableWithTotals( @tableHtml, @totalTableRow, 'confluence' );
    }
    else {
        $tableHtmlStr = "<b>Not available</b>";
    }

    my %values      = ();
    $values{projectName} = $projectName;

    # check all the files exist
    my $htmlFileName = "indel.html";
    my $htmlFilePath = File::Spec->catfile( $htmlPath, $htmlFileName );
    $values{title} = "$projectName - analysis: indels for all chromosomes";
    $values{path}  =
        "<a href=\"Home.html\">$projectName</a>"
      . "&gt;&nbsp;<a href=\"$htmlFileName\">indels for all chromosomes</a>";
    $values{table} = $tableHtmlStr;
    $values{page} = $htmlFileName;
    my $html = fromHtmlTemplate2HtmlFile( %values, 'body', $htmlFilePath );
    htmlFileMenuRegister( $htmlPath , $htmlFileName, "Ca1.7: Indels");
}
