package ScanFunctions;
use strict;
use ScanObject;

#This class builds a data structure that holds scanace information
#and performs calculations on the data
sub new {
    my $proto = shift;
    my $class = ref($proto) || $proto;
    my $self  = {};
    
    $self->{GENOME} = (); #Holds the names of all orfs in the genome of interest
    $self->{GENES} = (); #Holds the names of all the genes found found with scanace
    $self->{GENOMESITES} = (); #Count of all of the sites in the genome
    $self->{MOTIF} = (); #Holds the names of the genes found in the motif
    $self->{TESTSET} = (); #Holds the names of all the sequences that were used to find the motifs
    $self->{INTERSECTION} = (); #The genes found in both the motif and genome hits
    $self->{SITES} = (); #A hash that contains all of the genes in the genome
    
    $self->{GENOMEHIT} = 0; #Number of genes found with scanace
    $self->{GENOMEHITS} = 0;  #Number of total hits found in the genome with scanace
    $self->{GENOMEFREQ} = 0; #The frequency of the the motif in the genome
    $self->{GENOMEAVE} = 0; #The average number of sites per sequence
    $self->{GENOMESD} = 0; #The standard devation of the number of sites per sequence
    
    $self->{TESTSETHIT} = 0; #Number of genes found in the test set with scanace
    $self->{TESTSETHITS} = 0; #Number of total hits found in the test set
    $self->{TESTSETSITES} = 0; #Count of all of the possible binding sites in the test set
    $self->{TESTSETFREQ} = 0; #The frequency of the motif in the test set
    $self->{TESTSETAVE} = 0; #The average number of sites per sequence in the test set
    $self->{TESTSETSD} = 0; #The standard devation of the number of sites per sequence
    
    $self->{RATIO} = 0; #The ratio of testset motif frequency to genome motif frequency
    $self->{SUPPORT} = 0; #Log likelihood score
    $self->{GROUP} = 0; #The group specificity score (testsethit, $#genes + 1, $#intersection + 1, $#genome +1);
    $self->{SITE} = 0; #The site specificity score ($#genes+1, testsetsites, $#intersection+1, genomesites)
    
    $self->{WIDTH} = ''; #The width of the motif
    bless($self, $class);
    return $self;
};

sub DESTROY {
    my $self = shift;
    return;
};

sub load_genes{
    my $self = shift;
    $self->{MOTIF} = shift;
    return;
};

sub get_testhits{
    my $self = shift;
    if (defined $self->{TESTSETHITS}){
        return $self->{TESTSETHITS};
    }
    else {return 0}
};

sub set_testhits{
    my $self = shift;
    $self->{TESTSETHITS} = shift;
};

sub get_testave{
    my $self = shift;
    if (defined $self->{TESTSETAVE}){
        my $rounded = sprintf("%.2g", $self->{TESTSETAVE});  
        return $rounded
    }else{return 0};
};

sub set_testave{
    my $self = shift;
    $self->{TESTSETAVE} = shift;  
    return;
};

sub get_testsd{
    my $self = shift;
    if (defined $self->{TESTSETSD}){
        my $rounded = sprintf("%.2g", $self->{TESTSETSD});  
        return $rounded;
    }
    else {return 0};
};

sub set_testsd{
    my $self = shift;
    $self->{TESTSETSD} = shift;  
    return;
};

sub get_genomehits{
    my $self = shift;
    if (defined $self->{GENOMEHITS}){
        return $self->{GENOMEHITS};
    }
    else {return 0};
};

sub set_genomehits{
    my $self = shift;
    $self->{GENOMEHITS} = shift;
    return;
};

sub get_genomeave{
    my $self = shift;
    if (defined $self->{GENOMEAVE}){
        my $rounded = sprintf("%.2g", $self->{GENOMEAVE});  
        return $rounded;
    }
    else {return 0};
};

sub set_genomeave{
    my $self = shift;
    $self->{GENOMEAVE} = shift;  
    return;
};

sub get_genomesd{
    my $self = shift;
    if (defined $self->{GENOMESD}){
        my $rounded = sprintf("%.2g", $self->{GENOMESD});  
        return $rounded;
    }
    else {return 0};
};

sub set_genomesd{
    my $self = shift;
    $self->{GENOMESD} = shift;  
    return;
};

sub get_ratio{
    my $self = shift;
    if (defined $self->{RATIO}){
        my $rounded = sprintf("%.2g", $self->{RATIO});  
        return $rounded;
    }
    else {return 0};
};

sub set_ratio{
    my $self = shift;
    $self->{RATIO} = shift;
    return;
};

sub get_site{
    my $self = shift;
    if (defined $self->{SITE}){
        my $rounded = sprintf("%.2g", $self->{SITE});  
        return $rounded;
    }
    else {return 0};
};

sub set_site{
    my $self = shift;
    $self->{SITE} = shift;  
    return;
};

sub get_group{
    my $self = shift;
    if (defined $self->{GROUP}){
        my $rounded = sprintf("%.2g", $self->{GROUP});  
        return $rounded;
    }
    else {return 0};
};

sub set_group{
    my $self = shift;
    $self->{GROUP} = shift;  
    return;
};

sub get_support{
    my $self = shift;
    if (defined $self->{SUPPORT}){
        return $self->{SUPPORT};
    }
    else {return 0};
};

sub read_sites {
    my ($self, $file) = @_;
    my @array = ();
    open(IN, "< $file");
    $/ = 'Sites:';
    while (<IN>){push @array, [split(/\n/, $_)]};
    for my $i (1 .. $#{$array[1]}){
        my @list = split(/\t/, $array[1][$i]);
        if ($list[1] !~ /\d+/){next};
        if ($list[2] !~ /\d+/){next};
        if ($list[3] !~ /[AGTCN]+/){next};
        my $element = ScanObject->new;
        $element->put_element( $list[1], $list[2], $list[3], $list[4] );
        push @{$self->{SITES}{$list[0]}{ELEMENT}}, $element;
        $self->{SITES}{$list[0]}{COUNT} = $#{$self->{SITES}{$list[0]}{ELEMENT}} + 1;
        $self->{WIDTH} = length($list[3]);
    };
    close IN;
    $/ = "\n";
    return;
};

sub get_sites {
    my $self = shift;
    my %empty = ();
    if (defined $self->{SITES}){
        return %{$self->{SITES}};
    }
    else {return %empty};
};

#This function takes the test set sequence file and counts the number of possible sites
#while also recording all of the genes in the testset
sub get_testset{
    my ($self, $file) = @_;
    open (FILE, "< $file");
    my $line = undef;
    my $total = 0;
    my %holder = ();
    while (<FILE>){
        $line = $_;
        if ($line =~ /\>(NM\_\d{6})\_\d+/){$holder{$1}++; next};
        if ($line =~ /\>(\S+)/){ $holder{$1}++; next };
        chomp $line;
        $line =~ s/([AGTC])(N)/$1\t$2/g;
        $line =~ s/(N)([AGTC])/$1\t$2/g;
        my @fragments = split(/\t/, $line);
        for my $i (0 .. $#fragments){
            if (length($fragments[$i]) < $self->{WIDTH}){next};
            if ($fragments[$i] =~ /[AGTC]+/){
                $total = $total + length($fragments[$i]) - $self->{WIDTH} + 1;
            }
            else {next};
        };
    };
    close FILE;
    @{$self->{TESTSET}} = keys %holder;
    $self->{TESTSETSITES} = $total;
    return;
};

#This function counts all the possible motif binding sites within the genome while also recording
#all of the genes within the genome
sub get_genome{
    my ($self, $file) = @_;
    open (FILE, "< $file");
    my $line = undef;
    my $total = 0;
    my %holder = ();
    while (<FILE>){
        $line = $_;
        if ($line =~ /\>(NM\_\d{6})\_\d+/){$holder{$1}++; next};
        if ($line =~ /\>(\S+)/){ $holder{$1}++; next };
        chomp $line;
        $line =~ s/([AGTC])(N)/$1\t$2/g;
        $line =~ s/(N)([AGTC])/$1\t$2/g;
        my @fragments = split(/\t/, $line);
        for my $i (0 .. $#fragments){
            if (length($fragments[$i]) < $self->{WIDTH}){next};
            if ($fragments[$i] =~ /[AGTC]+/){
                $total = $total + length($fragments[$i]) - $self->{WIDTH} + 1;
            }
            else {next};
        };
    };
    close FILE;
    @{$self->{GENOME}} = keys %holder;
    $self->{GENOMESITES} = $total;
    return;
};

sub get_statistics {
    my $self = shift;
    my %holder = ();
    for my $i (0 .. $#{$self->{TESTSET}}){
        if (defined $self->{SITES}{$self->{TESTSET}[$i]}{COUNT}){
            if ($self->{SITES}{$self->{TESTSET}[$i]}{COUNT} > 0){
                $self->{TESTSETHIT} = $self->{TESTSETHIT} + 1;
            };
            $self->{TESTSETHITS} = $self->{TESTSETHITS} + $self->{SITES}{$self->{TESTSET}[$i]}{COUNT};
        };
    };
    $self->{TESTSETAVE} = $self->{TESTSETHITS}/($#{$self->{TESTSET}} + 1);
    $self->{TESTSETFREQ} = $self->{TESTSETHITS} / $self->{TESTSETSITES};
    my $testsd = 0;
    for my $i (0 .. $#{$self->{TESTSET}}){
        if (defined $self->{SITES}{$self->{TESTSET}[$i]}{COUNT}){
            $testsd = $testsd + (($self->{SITES}{$self->{TESTSET}[$i]}{COUNT} - $self->{TESTSETAVE}) ** 2);
        }
        else {$testsd =  $testsd + ($self->{TESTSETAVE} ** 2)};
    };
    $self->{TESTSETSD} = sqrt($testsd/$#{$self->{TESTSET}});
    for my $i (0 .. $#{$self->{GENOME}}){
        if (defined $self->{SITES}{$self->{GENOME}[$i]}{COUNT}){
            if ($self->{SITES}{$self->{GENOME}[$i]}{COUNT} > 0){
                $holder{$self->{GENOME}[$i]}++;
                $self->{GENOMEHIT} = $self->{GENOMEHIT} + 1;
            };
            $self->{GENOMEHITS} = $self->{GENOMEHITS} + $self->{SITES}{$self->{GENOME}[$i]}{COUNT};
        };
    };
    my @hits = keys %holder;
    $self->{GENES} = \@hits;
    $self->{GENOMEAVE} = $self->{GENOMEHITS}/($#{$self->{GENOME}} + 1);
    my $genomesd = 0;
    for my $i (0 .. $#{$self->{GENOME}}){
        if (defined $self->{SITES}{$self->{GENOME}[$i]}{COUNT}){
            $genomesd = $genomesd + (($self->{SITES}{$self->{GENOME}[$i]}{COUNT} - $self->{GENOMEAVE}) ** 2);
        }
        else { $genomesd = $genomesd + ($self->{GENOMEAVE} ** 2)};
    };
    if ($#{$self->{GENOME}} <= 0){$self->{GENOMESD} = 0} 
    else{ $self->{GENOMESD} = sqrt($genomesd/$#{$self->{GENOME}})};
    if ($self->{GENOMESITES} <= 0){$self->{GENOMEFREQ} = 0}
    else{ $self->{GENOMEFREQ} = $self->{GENOMEHITS} / $self->{GENOMESITES}};
    if ($self->{GENOMEFREQ} <= 0){$self->{RATIO} = 0}
    else{ $self->{RATIO} = $self->{TESTSETFREQ} / $self->{GENOMEFREQ}};
    if ($self->{TESTSETFREQ} > 0 and $self->{GENOMEFREQ} > 0){
        $self->{SUPPORT} = ($self->{TESTSETHITS} * log($self->{TESTSETFREQ})) - ($self->{GENOMEHITS} * log($self->{GENOMEFREQ}))
    }
    else{$self->{SUPPORT} = 0};
    my %union = ();
    my %intersect = ();
    foreach my $e (@{$self->{GENES}}, @{$self->{MOTIF}}){$union{$e}++ && $intersect{$e}++};
    my @intersect = keys %intersect;
    $self->{INTERSECTION} = \@intersect;
    $self->{GROUP} = hypergeometric(($#{$self->{INTERSECTION}} + 1), ($#{$self->{TESTSET}} + 1), ($#{$self->{GENES}} + 1), ($#{$self->{GENOME}}+1));
    $self->{SITE} = hypergeometric(($#{$self->{INTERSECTION}} +1 ), ($#{$self->{GENES}} + 1), $self->{TESTSETSITES}, $self->{GENOMESITES});
    return;
};

sub clean {
    my $self = shift;
    my %holder = ();
    for my $i (0 .. $#{$self->{TESTSET}}){
        if (exists $self->{SITES}{$self->{TESTSET}[$i]}){
            $holder{$self->{TESTSET}[$i]}{COUNT} = $self->{SITES}{$self->{TESTSET}[$i]}{COUNT};
            for my $j (0 .. $#{$self->{SITES}{$self->{TESTSET}[$i]}{ELEMENT}}){
                my $temp = new ScanObject;
                $self->{SITES}{$self->{TESTSET}[$i]}{ELEMENT}[$j]->copy($temp);
                push @{$holder{$self->{TESTSET}[$i]}{ELEMENT}}, $temp;
            };
        };
    };
    $self->{SITES} = \%holder;
    $self->{GENOME} = ();
    return;
};

sub lngamma {
    my $c = shift;
    my @cof = (76.18009172947146,
             -86.50532032941677,
             24.01409824083091,
             -1.231739572450155,
             0.1208650973866179e-2,
             -0.5395239384953e-5);
    my $xx = $c;
    my $yy = $c;
    my $tmp = $xx + 5.5 - ($xx + 0.5) * log($xx + 5.5);
    my $ser = 1.000000000190015;
    for(my $j = 0; $j<=5; $j++){
        $ser += ($cof[$j] / ++$yy)};
    return(log(2.5066282746310005*$ser/$xx)-$tmp);
}

sub lnbico {
    my $n = shift;
    my $k = shift;
    return (lngamma($n+1) - lngamma($k+1) - lngamma($n-$k+1));
};

sub hypergeometric {
    my $x = shift; #Intersection of motif sequences and scanace hits
    my $s1 = shift; #Total sequences used to make the motif
    my $s2 = shift; #Total sequences identified by scanace as hits
    my $N = shift; #Total sequences in the genome
    my $lnhyper = lnbico($s1, $x) + lnbico(($N-$s1) , ($s2-$x)) - lnbico($N, $s2);
    return exp($lnhyper);
};

sub analysis {
    my ( $self, $scanfile, $sequences, $genome, $genearray ) = @_;
    $self->{MOTIF} = $genearray;
    $self->read_sites($scanfile);
    $self->get_testset($sequences);
    $self->get_genome($genome);
    $self->get_statistics();
    $self->clean();
    return;
};

1;
