#!/usr/bin/env perl

=begin
analyses patient sample vcf files for variant frequencies across multiple runs.
Produces 1 table + 1 graph for each sample in each of n = 1, 2, 3, etc directories.
n values correspond to frequency each variant found across multiple runs, so in
n = 1, all html tables list variants found only once for that sample across all
runs analysed. In n = 2, tables list variants found in 2 out of total number of
runs, etc. Each variant can appear >= once per vcf file (html table).
=cut

use lib '/home/raj/perl5/lib/perl5';

use Modern::Perl '2012';
use Data::Dumper::Concise;
use File::Path qw(make_path remove_tree); # IO::All PITA for dir manipulation
use Path::Tiny;
use FindBin; # warn $FindBin::Bin;
use IO::All;
use CGI;

use lib $FindBin::RealBin . '/../../lib';
use NGS::ChartMaker;
use DBIx::Simple;

my $src_dir = $FindBin::Bin . '/src/analyse';
my @dirs = io($src_dir)->all; # print "$_ is a " . $_->type . "\n" for @dirs; exit;
    
my %samples;
# collect filenames from each dir, add to %samples hash:
for my $d (@dirs) { # warn Dumper $_;
    next unless $d->type eq 'dir'; # skip eg png's
	my @files = io($d)->all; # Dumper \@files;  
	for my $f (@files) { # warn $f->filename;
		my $filename = $f->filename;
		next unless $filename =~ /vcf\Z/; # warn $f->filename;
		push @{ $samples{$filename} }, $f; # io object
	}
} # warn Dumper \%samples; exit;

# delete and recreate analysis output dirs, one for each src dir:
for ( my $i = 1; $i <= @dirs; $i++ ) { # warn $i;
    my $dir = $FindBin::Bin . '/analysis/' . $i;
    remove_tree($dir) if -e $dir;  # not worth the hassle using IO::All   
    make_path($dir); # ditto
}

my $dbix = _build_dbix(); # warn Dumper $dbix;
my $EXCLUDES = _get_excludes(); # warn Dumper $EXCLUDES;

my $q = CGI->new;
my %h = (); # start-position indexer (global, used in parse(), reset in FILE loop)

my %chart = (
    x_title => 'Depth',
    y_title => 'VF',
);
my %all_data;

# params to collect from vcf file:
my @params = qw( position vf alt ref gene filter depth runId);

FILE: # check each file is represented in all dirs:
while ( my($sample, $o) = each %samples ) {
	unless ( scalar @$o == scalar @dirs ) { # one object per file per directory
		printf "%s seen only %s times\n", $sample, scalar @$o;
		next FILE;
	}
    
	%h = (); # reset
	# parse file content in each directory, %h holds array of data
	for my $f(@$o) { # warn $f->name;
		parse($f);
	} # warn Dumper \%h;
	
	my %data; # array of all results (href) for this sample:
	while ( my($position, $d) = each %h ) { # warn Dumper $d;
		my ( $n, $results ) = ( $d->{count}, $d->{data} );

		for my $result(@$results) { 
			# add position to $result href:
			$result->{position} = $position;
			push @{ $data{$n} }, $result;
		}
	}

    $sample =~ s/\.vcf//; # warn $sample;
    
    # now have hash of all results from all dirs for this sample, keys = $n (1,2,3,etc)
    while ( my($n, $result) = each %data ) {
        my @rows = (); # reset
        $chart{data} = {}; # reset
        # sort each row by variation frequency (highest -> lowest):
        for my $row ( sort by_var_freq @$result ) {
            my $position = $row->{position};
            
            push @rows, $q->Tr( $q->td( [ @{$row}{@params} ] ) );
            
            push @{ $chart{data}{x} }, $row->{depth};
			push @{ $all_data{$n}{x} }, $row->{depth};
			
            push @{ $chart{data}{y} }, $row->{vf};
			push @{ $all_data{$n}{y} }, $row->{vf};
        } # warn Dumper \@rows; exit;

        my $total_variants = do {
            my %h; $h{ $_->{position} }++ for @$result;
            scalar keys %h;
        };
        my $content = do {
            my %d = (
                sample => $sample,
                total  => $total_variants,
                rows   => \@rows,
            );
            _template(\%d);
        };

        my $file = sprintf './analysis/%s/%s.html', $n, $sample;
        $content > io($file);        
        
        $chart{chart_title} = $sample; # warn Dumper scalar @{ $chart{data}{x} };
        $chart{data_title}  = sprintf '%s/%s', $n, $sample;
        NGS::ChartMaker->new(%chart)->make_chart();
    }
}

# output combined graphs, one each for n=1, n=2:
for (1,2) {
	my %h = (
		x_title => 'Depth',
		y_title => 'VF',
		chart_title => "All data $_",
		data_title => 'all_data_'.$_,
		data => $all_data{$_},
	); # warn Dumper scalar @{ $h{data}{y} };
	NGS::ChartMaker->new(%h)->make_chart();
}

# expand x-axis scale (0-4000) on n=1:
{
	my %h = (
		x_title => 'Depth',
		y_title => 'VF',
		max_x_val => 4000,
		chart_title => 'All data #1 expanded',
		data_title => 'all_data_1_expanded',
		data => $all_data{1},
	);
	NGS::ChartMaker->new(%h)->make_chart();
}

exit 0;

#===============================================================================

sub parse {
	my $f = shift; # io object
    
	my @folders = $f->splitdir; # warn $folders[-2]; # array of dirs in full path
	my @lines   = io($f->name)->slurp; # warn Dumper \@lines;
    
	my %local = (); # indexer localised to individual file

	LINE: for (@lines) {
		next LINE if /^#/; # skip comments
		my @fields = split "\t", $_;
		
		my $start_pos = $fields[1]; # warn $start_pos;
			# warn $start_pos if $EXCLUDES->{$start_pos};
		next LINE if (! $start_pos) || $EXCLUDES->{$start_pos};  
			# warn $start_pos if $EXCLUDES->{$start_pos};

        # initialise intra-file start position indexer:
		$local{$start_pos} ||= 1; # only need true/false, not item frequency
		
        { # add data as array(ref) - can have multiple start pos entries per file:
            my $data = _parse_fields(\@fields);
            # add run directory name:
            $data->{runId} = $folders[-2]; # warn $data->{runId}; # 2nd to last
            push @{ $h{$start_pos}{data} }, $data;
        }
	}
	# auto-increment %h start position key for each start position key in %i:
	for my $start_position (keys %local) { $h{$start_position}{count}++ }
}

sub by_var_freq {
	return $b->{vf} <=> $a->{vf} || $a->{depth} <=> $b->{depth};
}

sub _template {
    my $args = shift; # sample name (str), total (int), rows (arrayref)
    
    my ($sample, $total, $rows) = @{$args}{qw/sample total rows/};
    
    return sprintf <<"HTML", join "\n", @$rows;
    <!DOCTYPE html>
    <html lang="en">
      <head>
        <style>
            html { font-family: verdana, sans-serif; }
            h2 { color: #800000; }
            h3 { }
            table { border-collapse: collapse; margin-left: 1em }
            td, th { border: solid 1px #c0c0c0; padding: 2px; font-size: 9pt; }
            #head { background-color: ActiveCaption; color: CaptionText; }
        </style>
      </head>
      <body>
        <h2>${sample}.vcf</h2>
        <h3>Total variants: $total</h3>
        <table>
          <tr id="head">
            <th>position</th>
            <th>vf</th>
            <th>alt</th>
            <th>ref</th>
            <th>gene</th>
            <th>filter</th>
            <th>depth</th>
            <th>run Id</th>
          <tr>
          <tbody>%s</tbody>
        </table>
      </body>
    </html>
HTML
}

sub _build_dbix {
    my $db = path($FindBin::RealBin, '..', '..', 'ngs.sqlite')->realpath;
	my $dsn = "dbi:SQLite:dbname=$db"; # warn 'DB:'. $db;

	my $dbix = DBIx::Simple->connect($dsn); # warn Dumper $dbix;
    return $dbix;
}

sub _get_excludes {
	$dbix->select('locations', ['position', 1], {action => 'exclude'})->map;
}

sub _parse_fields {
    my $data = shift; # arrayref of tab-separated vcf row entries

	my $ref     = $data->[3]; # warn $ref;
	my $alt     = $data->[4]; # warn $alt;
	my $filter  = $data->[6]; # warn $filter;
    my $info    = $data->[7]; # warn $info
	my $results = $data->[9]; # warn $results;
    
	my $vf = (split ':', $results)[3]; # warn $vf;
	my ($depth) = $info =~ /DP=(\d+)/; # warn $depth;
    my ($genes) = $info =~ /GI=(.*?)\;/; # warn $genes; # non-greedy match

    my $gene = do { # remove duplicate genes: eg TP53,TP53,TP53
        no warnings 'uninitialized'; # eg if regex failed due to malformed str
        my %seen = map { $_ => 1 } split ',', $genes;    
        join ',', keys %seen;
    } || '[UNDEF]';
    my %h = (
		filter => $filter,
		depth  => $depth,
        gene   => $gene,
		ref    => $ref,
		alt    => $alt,
		vf     => $vf,
	);
    return \%h;
}