package NGS::VEP;

# wrapper around Ensembl Variant Effect Predictor stand-alone script;
# http://www.ensembl.org/info/docs/variation/vep/vep_script.html

use Moo;
use IO::All;
use IO::File;
use Text::CSV;
use Data::Dumper;
use FindBin qw($Bin); # warn $Bin;

use MooX::Types::MooseLike::Base qw(:all);
use NGS::MooX::Types qw(HashReference);
use NGS::DB;

# required to be passed on object contruction:
has src_data  => ( is  => 'ro', isa => Str, required => 1 );
has form_opts => ( is  => 'ro', isa => HashReference, required => 1 );

# called internally by result():
has ref_table => ( is => 'lazy' ); # _build_ref_table()

my %text_csv_args = ( sep_char => "\t", binary => 1 ); # global args for Text::CSV

# location of variable_effect_predictor script:
my $vep_script = "$Bin/../script/variant_effect_predictor.pl";

#===============================================================================
sub result {
    my $self = shift; # warn Dumper $self->database;
    
    my $data = $self->_text_csv_slurp(); # warn Dumper $data; # arrayref

    # get reference table:
    my $ref_table = $self->ref_table; # warn Dumper $reference;
    
    # create vep file for variant_effect_predictor.pl
    my $vep_file = "$Bin/../results.vep";
    io($vep_file)->unlink if -e ($vep_file); # clear previous

    my @vep_fields = qw(chromosome start_point end_point allele strand exon);    
    my %vep_data; # will hold ref data for merging with VEP script output
    my $i; # will hold total count of 'accepted' rows

    ROW: for my $row (@$data) { # cols = ref_seq, base/variant, status, ....
        next ROW unless lc $row->[2] eq 'accepted'; # 3rd col is status field
        $i++; # increment 'accepted' row counter

        my $exon = $row->[0]; # warn $exon;

        # get data from reference table, or next record:
        my $ref = $ref_table->{$exon}; # warn Dumper $ref;
        if (! $ref) {
            warn "no ref table entry for $exon"; next ROW;
        }
        my $chromosome = $ref->{chromosome};
        
        # split 2nd col to get location & variant (eg 136:A/T, 26-27:TG/--)
        my ($location, $allele) = split ':', $row->[1]; # warn Dumper [$location, $allele];
        
        # define initial start & end point as ref table start position minus 1:
        my $start_point = my $end_point = $ref->{start_base} - 1;

        # if location > 1 base (eg 26-27, 134-137 , etc):
        if ( $location =~ /(\d+)-(\d+)/ ) { # warn $location; # eg 26-27
			my $span = $2 - $1; # eg 26-27 will give span = 1
            $start_point += $1; # start location (26 in above example)
            $end_point = $start_point + $span; 
        }
        else {
            $end_point = $start_point += $location;
        }
        
        { # manipulate exon name to provide unique ID for vep input
          # eg TET2 exon 4 [891], TET2 exon 4 [1031], etc:
            $exon =~ s/\s/%%/g; # vep can't deal with spaces in any col
            $exon .= '%%[' . $location . ']'; # exon + location gives unique ID
        }
        { # data structure for vep script input file: 
            my %h = (
                chromosome  => $chromosome,
                start_point => $start_point,
                end_point   => $end_point,
                allele      => $allele,
                strand      => '+', # required 5th col in VEP file
                exon        => $exon, # unique identifier eg "TET2 exon 4 [891]"
            ); # warn Dumper \%h;
            
            # append to vep script input file:
            my $line = join "\t", @h{@vep_fields};
            io($vep_file)->append($line . "\n");
            
            { # data for use with vep output (use chr:start_location as unique ID):
              # vep returns eg 2:12345 or 2:12345-12347 in 'Location' col
                my $base_num = ( $start_point == $end_point )
                    ? $start_point # eg 2:12345
                    : $start_point .'-'. $end_point; # eg 2:12345-12347    
                my $unique_id = join ':', $chromosome, $base_num;
                $vep_data{$unique_id} = \%h;
            }
        }
    } # warn Dumper \%vep_data;
   
    my $result = $self->_run_vep(\%vep_data); # hashref of data & orphaned rows    
    $result->{row_count} = $i; # add row count of accepted rows 
    return $result;
}

#-------------------------------------------------------------------------------
sub _run_vep {
    my ($self, $src_data) = @_; # HoH "chr:start" => (start_point, end_point, allele, etc)
    
	my @results; # to hold tab-delimited results array from vep output

    # load opts from vep.ini + user-selected options to augment:
    my @opts = ( "-config=$Bin/../script/vep.ini" );
    my $user_opts = $self->_get_form_options;
    push @opts, @$user_opts;
    
    # open pipe to vep script (input = $vep_file; output = stdout):
	open my $fh, "-|", $vep_script, @opts;
	while ( my $line = readline($fh) ) { # warn $line; 
		push @results, $line;
	} # warn Dumper \@results;
    
    my @input_fields = qw(exon chromosome start_point end_point allele);
    
    my @data; # to hold composite of vep input & vep output
    my %seen; # running tally of vep inputs with 1 or more results

    for (@results) {
        next if /^#/; # skip comment lines
        my @cols = split /\t/;
        
        my $chr_location = $cols[1]; # in same format as %$src_data key (eg 2:123456)
        
        if ( my $ref = $src_data->{$chr_location} ) {
            # capture vep cols: consequence, amino_acids, existing_variation, extra:
            my @output_data = @cols[6,10,12,13]; 
            push @data, [ @{$ref}{@input_fields}, @output_data ]; # combine data
            $seen{$chr_location}++; # vep input file has at least one result
        }
        else { warn $_ } # no match in src file
    }
    
    # has no vep result:
    my @orphans = sort by_chr_location
        map $src_data->{$_}, grep { ! $seen{$_} } keys %$src_data;
    
    return { data => \@data, orphans => \@orphans }
}

#-------------------------------------------------------------------------------
sub _get_form_options { # user submitted opts (eg sift, polyphen, etc)
    my $self = shift;

    my $params = $self->form_opts; # warn Dumper $form_opts;
    
    my @opts;
    
    # polyphen & sift take args (s, p or b)
    push @opts, '--polyphen=' . $params->{polyphen} if $params->{polyphen};
    push @opts, '--sift='     . $params->{sift} if $params->{sift};
    # check_existing, regulatory & coding_only take no args:
    for ( qw/check_existing coding_only regulatory/ ) {
        push @opts, "--$_" if $params->{$_};
    } # warn Dumper \@opts;
    
    return \@opts;
}

#-------------------------------------------------------------------------------
sub by_chr_location { # alphanumeric sort on chromosome 1-22, X, and start_point
    return ( # 1st sort on chromosome, then start_point:
        ( $a->{chromosome} !~ /^\d+/ || $b->{chromosome} !~ /^\d+/ )
            ? $a->{chromosome} cmp $b->{chromosome} # non-digit eg X
            : $a->{chromosome} <=> $b->{chromosome} # digit eg 1 - 22
    ) || $a->{start_point} <=> $b->{start_point};
}

#-------------------------------------------------------------------------------
# extract tab-delimited data from csv source file:
sub _text_csv_slurp {
    my $self = shift;
    
    my $data = $self->src_data; # string

    my $csv = Text::CSV->new(\%text_csv_args);
    
    my @results;
    for ( split "\n", $data ) { # warn $_;
        if ( $csv->parse($_) ) {
            my @row = $csv->fields();
            push @results, \@row;
        }
    }
    return \@results;
}

#-------------------------------------------------------------------------------
# get reference table as hashref of gene => { chromosome, start_base }:
sub _build_ref_table {
    my $self = shift;
    
    my $db = NGS::DB->new(); # warn $db->dbix;
    my $ref_tbl = $db->dbix->select('ref_seq', '*')
        ->map_hashes('gene'); # warn Dumper $ref_tbl;
    return $ref_tbl;
}
 
=begin _build_ref_table() from ref.csv
    my $ref = "$Bin/../ref.csv";
    my $io = new IO::File;
    open( $io, '<', $ref ) || die $!;

    my $csv = Text::CSV->new(\%text_csv_args);

    my %h;
    while ( my $row = $csv->getline( $io ) ) {
        my ($ref_seq, $chr, $start) = @$row;
        $h{$ref_seq} = {
            chromosome => $chr,
            start_base => $start,
        };
    } # warn Dumper \%h;
    
    return \%h;
=cut

1;