#!/usr/bin/env perl
# PROJECT: CASAVA
# MODULE:  $RCSfile: generateVariantsReports.pl,v $
# AUTHOR:  Chris Saunders
#
# 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 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 Casava::Common::Log;
use Casava::PostAlignment::Sequencing::Config qw(loadConfiguration
  %CONF_APP %CONF_PROJ);
use Casava::Common::IOLib qw(file2tableEx);
use Casava::Common::Template;
use Casava::PostAlignment::Sequencing::Reports qw(getColIndex htmlFileMenuRegister makeChart updateTotalRow);


my $argvStr     = join ' ', @ARGV;
my $cmdline     = "$0 $argvStr";
my $scriptName = (File::Spec->splitpath($0))[2];

my $usage =
  "$scriptName [options]\n" .
  "\t--projectDir - Project directory\n"
;
my $projectDir = "";
my $help       = 0;
my $sampleRateSetSuffix = '';
my $result = GetOptions(
                         'projectDir|p=s' => \$projectDir,
                         'sampleRateSetSuffix=s' => \$sampleRateSetSuffix,
                         'help|h'         => \$help
);

if ((not $result) or $help) {
    print "$usage\n\n";
    exit(0);
}

errorExit "ERROR: Incorrect number of arguments\n$usage" if ( $projectDir eq "" );

errorExit "ERROR: Unrecognized arguments: " . join(" ",@ARGV) . "\n" if(@ARGV);


{
    my $hostname = hostname();
    printLog( "Running $hostname:[$cmdline]\n", 0);
}

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

my $isNoIndels      = ((defined $CONF_PROJ{variantsNoIndels}) and $CONF_PROJ{variantsNoIndels});

my $htmlPath  = File::Spec->catdir( $projectDir, "html" );
my $statsPath = File::Spec->catdir( $projectDir, "stats" );

{
    my $dataDir = '/home/psgendb/local/pkg/CASAVA_v1.8.2-build/share/CASAVA-1.8.2';
    my $templPath = File::Spec->catdir( $dataDir, "templates" );
    configure($templPath);
}



sub makePage($$$$) {
    my ($htmlFileName,$description,$tableHtmlStrRef,$label) = @_;

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

    my $htmlFilePath = File::Spec->catfile( $htmlPath, $htmlFileName );
    $values{title} = "$projectName - analysis: $description";
    $values{path}  =
    	"<a href=\"Home.html\">$projectName</a>"
      . "&gt;&nbsp;<a href=\"$htmlFileName\">$description</a>";
    $values{table} = $$tableHtmlStrRef;
    $values{page} = $htmlFileName;
    fromHtmlTemplate2HtmlFile( %values, 'body', $htmlFilePath );
    htmlFileMenuRegister( $htmlPath , $htmlFileName, $label);
}



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



sub formatCovRow {
    my ($rowRef,$colRef) = @_;

    my $sites = $$rowRef[ $$colRef{sites} ];
    my $known_sites = $$rowRef[ $$colRef{known_sites} ];
    my $known_sites_mapped = $$rowRef[ $$colRef{known_sites_mapped} ];
    my $known_sites_coverage = $$rowRef[ $$colRef{known_sites_coverage} ];

    return ($$rowRef[ $$colRef{seq_name} ],
            $sites,
            $known_sites,
            $known_sites_coverage,
            sprintRatio($known_sites_coverage,$known_sites),
            sprintRatio($known_sites_mapped,$known_sites));
}



# generate html report for coverage:
#
if ( 1 ) {
    my @tableHeader = (
                        'name',
                        'sites',
                        'known sites',
                        'bases mapped at known sites',
                        'mean depth at known sites',
                        'fraction of known sites mapped'
                      );
    my @tableHtml = ( \@tableHeader );

    my $tableHtmlStr = '';
    my $fileIn   = File::Spec->catfile( $statsPath, "coverage.summary.txt" );

    if ( -e $fileIn ) {
        my @Xaxis                  = ();
        my @meanDepthData          = ();
        my @fractionMappedData     = ();

        my @tableIn = ();
        my %fileMD = ();
        file2tableEx( @tableIn, %fileMD,  $fileIn );

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

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

            push @Xaxis,              $rowRef->[ $colIndex{seq_name} ];
            push @meanDepthData,      $tableRow[4];
            push @fractionMappedData, $tableRow[5];
        }

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

        # copying high depth exclusion routine from CASAVA 1.7:
        #
        #chrom M value is typically so high that all others will stop being
        #visible if we plot it.
        my @Xaxis_filter;
        my @meanDepthData_filter;
        {
            my $maxPlotDepthFactor = 4;
            my $maxPlotMeanDepth = $totalTableRow[4]*$maxPlotDepthFactor;
            for(my $i=0;$i<scalar(@Xaxis);++$i){
                next if($meanDepthData[$i] > $maxPlotMeanDepth);
                push @Xaxis_filter,$Xaxis[$i];
                push @meanDepthData_filter,$meanDepthData[$i];
            }
        }

        {
            my $plotTitle = "mean depth at known sites";
            my $pngFile = "coverage.mean.png";
            my $xAxisLabel = "sequence name";
            my $yAxisLabel = "mean depth at known sites";

            makeChart($htmlPath,
                      $pngFile,
                      \$tableHtmlStr,
                      $plotTitle,
                      $xAxisLabel,
                      $yAxisLabel,
                      \@Xaxis_filter,
                      [ [ \@meanDepthData_filter,   0x80ff80, "depth" ] ]);
        }

        {
            my $plotTitle = "fraction of known sites mapped";
            my $pngFile = "coverage.mapped.png";
            my $xAxisLabel = "sequence name";
            my $yAxisLabel = "fraction of known sites mapped";
            makeChart($htmlPath,
                      $pngFile,
                      \$tableHtmlStr,
                      $plotTitle,
                      $xAxisLabel,
                      $yAxisLabel,
                      \@Xaxis,
                      [ [ \@fractionMappedData,   0x80ff80, "mapped" ] ]);
        }

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

    my $htmlFileName = "coverage.html";
    my $description = "coverage for all reference sequences";
    my $label = "Coverage";
    makePage($htmlFileName,$description,\$tableHtmlStr,$label);
}



sub formatSnpRow {
    my ($rowRef,$colRef) = @_;

    my $hom = $$rowRef[ $$colRef{hom} ];
    my $het_ref = $$rowRef[ $$colRef{het} ];
    my $het_nonref = $$rowRef[ $$colRef{het_nonref} ];
    my $het = $het_ref + $het_nonref;
    my $total = $hom + $het;
    my $hetHomRatio = sprintRatio($het,$hom);

    return ($$rowRef[ $$colRef{seq_name} ],
            $het_ref,
            $het_nonref,
            $het,
            $hom,
            $hetHomRatio,
            $total);
}


# generate html report for snps:
#
if ( 1 ) {
    my @tableHeader = (
                        'name',
                        'ref het',
                        'nonref het',
                        'total het',
                        'total hom',
                        'het/hom ratio',
                        'total snps'
                      );
    my @tableHtml = ( \@tableHeader );

    my $tableHtmlStr = '';
    my $fileIn   = File::Spec->catfile( $statsPath, "snps.summary.txt" );

    if ( -e $fileIn ) {
        my @Xaxis                  = ();
        my @heterozygousData       = ();
        my @homozygousData         = ();

        my @tableIn = ();
        my %fileMD = ();
        file2tableEx( @tableIn, %fileMD,  $fileIn );

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

        my $minQsnp = $fileMD{SUMMARY_MIN_QSNP};

        my @totalRow = ();
        for my $rowRef (@tableIn) {
            updateTotalRow(\@totalRow,$rowRef,$colIndex{seq_name});

            push @Xaxis,                  $rowRef->[ $colIndex{seq_name} ];
            push @heterozygousData,       $rowRef->[ $colIndex{het} ] + $rowRef->[ $colIndex{het_nonref} ];
            push @homozygousData,         $rowRef->[ $colIndex{hom} ];

            my @tableRow = formatSnpRow($rowRef,\%colIndex);
            push @tableHtml, \@tableRow;
        }

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

        {
            my $plotTitle = "called snp type summary";
            $plotTitle .= " for Q(snp) >= $minQsnp" if(defined $minQsnp);
            my $pngFile = "snps.png";
            my $xAxisLabel = "sequence name";
            my $yAxisLabel = "number of snps";

            makeChart($htmlPath,
                      $pngFile,
                      \$tableHtmlStr,
                      $plotTitle,
                      $xAxisLabel,
                      $yAxisLabel,
                      \@Xaxis,
                      [ [ \@homozygousData,   0x80ff80, "homozygous" ] ,
                        [ \@heterozygousData, 0xff8080, "heterozygous" ] ]);
        }

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

    my $htmlFileName = "snps.html";
    my $description = "snps for all reference sequences";
    my $label = "SNPs";
    makePage($htmlFileName,$description,\$tableHtmlStr,$label);
}



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

    my $hom = $$rowRef[ $$colIndexRef{hom} ];
    my $het = $$rowRef[ $$colIndexRef{het} ];
    my $insertion = $$rowRef[ $$colIndexRef{insertion} ];
    my $deletion = $$rowRef[ $$colIndexRef{deletion} ];
    my $breakpoint = $$rowRef[ $$colIndexRef{breakpoint} ];
    my $totalIndel = $hom + $het;
    my $hetHomRatio = sprintRatio($het,$hom);

    return ($$rowRef[ $$colIndexRef{seq_name} ],
            $insertion,
            $deletion,
            $breakpoint,
            $hom,
            $het,
            $hetHomRatio,
            $totalIndel);
}



# generate separate html report for indels:
#
if (not $isNoIndels) {
    my @tableHeader = (
                        'name',
                        'insertion',
                        'deletion',
                        'breakpoint',
                        'homozygous',
                        'heterozygous',
                        'het/hom ratio',
                        'total indels'
    );
    my @tableHtml     = ( \@tableHeader );

    my $tableHtmlStr = '';
    my $fileIn   = File::Spec->catfile( $statsPath, "indels.summary.txt" );

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

        my @tableIn  = ();
        my %fileMD = ();
        file2tableEx( @tableIn, %fileMD,  $fileIn );

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

        my $minQindel = $fileMD{SUMMARY_MIN_QINDEL};

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

            push @Xaxis,                  $rowRef->[ $colIndex{seq_name} ];
            push @heterozygousData,       $rowRef->[ $colIndex{het} ];
            push @homozygousData,         $rowRef->[ $colIndex{hom} ];
            push @insertionData,          $rowRef->[ $colIndex{insertion} ];
            push @deletionData,           $rowRef->[ $colIndex{deletion} ];
            push @breakpointData,         $rowRef->[ $colIndex{breakpoint} ];

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

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

        my $xAxisLabel = "sequence name";
        my $yAxisLabel = "number of indels";

        {
            my $plotTitle = "called indel category summary";
            $plotTitle .= " for Q(indel) >= $minQindel" if(defined $minQindel);
            my $pngFile = "indels.cat.png";

            makeChart($htmlPath,
                      $pngFile,
                      \$tableHtmlStr,
                      $plotTitle,
                      $xAxisLabel,
                      $yAxisLabel,
                      \@Xaxis,
                      [ [ \@insertionData,   0x80ff80, "insertion" ],
                        [ \@deletionData, 0xff8080, "deletion" ],
                        [ \@breakpointData, 0x8080ff, "breakpoint" ] ]);
        }

        {
            my $plotTitle = "called indel type summary";
            $plotTitle .= " for Q(indel) >= $minQindel" if(defined $minQindel);
            my $pngFile = "indels.type.png";

            makeChart($htmlPath,
                      $pngFile,
                      \$tableHtmlStr,
                      $plotTitle,
                      $xAxisLabel,
                      $yAxisLabel,
                      \@Xaxis,
                      [ [ \@homozygousData,   0x80ff80, "homozygous" ] ,
                        [ \@heterozygousData, 0xff8080, "heterozygous" ] ]);
        }

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

    my $htmlFileName = "indels.html";
    my $description = "indels for all reference sequences";
    my $label = "Indels";
    makePage($htmlFileName,$description,\$tableHtmlStr,$label);
}

printLog( "$0 Finished\n", 0);

1;
