package SequenceFunctions;
use strict;
use Storable;
use MotifFunctions;
use XML::Writer;
use IO::File;

sub new {
    my $proto = shift;
    my $class = ref($proto) || $proto;
    my $self  = {};
    $self->{SETS}   = (); #Contains an hash of gene clusters pointing to arrays of gene names
    $self->{SEARCHSET} = (); #Contains a list of genes used to search for ScanACE sites
    $self->{GENOME} = (); #This data structure contains Annotation about genes
    $self->{CHROMOSOMES} =(); #This is an array that holds the names of all chromosomes used
    $self->{CONSERVED} = (); #This is where the locations of the conserved regions are stored
    $self->{NONCONSERVED} = (); #Regions of nonconserved DNA
    $self->{INTRONIC} = {}; #Intronic regions of dna
    $self->{EXONIC} = {}; #exonic regions of dna
    $self->{MOTIFS} = ();
    $self->{FACTORS} = ();
    $self->{FREQUENCY} = (); #holder for markov model
    $self->{SEQUENCES} = ();
    $self->{GENES} = 0;
    $self->{DIRECTORIES} = { GENOME => '',
                             WORK => '',
                             BLAST => '',
                             TEMP => ''
                           };
    $self->{INFO} = {UPSTREAM => '',
                     DOWNSTREAM => '',
                     CONSERVED => '0',
                     MASKED => '1',
                     MINIMUM => '0'
                    };
    bless($self, $class);
    return $self;
};

sub DESTROY {
    my $self = shift;
};

sub set_info {
    my $self = shift;
    $self->{INFO}{UPSTREAM} = shift;
    $self->{INFO}{DOWNSTREAM} =shift;
    $self->{INFO}{CONSERVED} = shift;
    $self->{INFO}{MASKED} = shift;
    $self->{INFO}{MINIMUM} = shift;
};

sub set_all_directories{
    my $self = shift;
    $self->{DIRECTORIES}{GENOME} = shift;
    $self->{DIRECTORIES}{WORK} = shift;
    $self->{DIRECTORIES}{TEMP} = shift;
    $self->{DIRECTORIES}{BLAST} = shift;
    print"Genome directory is located\:\n$self->{DIRECTORIES}{GENOME}\n";
    print"Working directory is located\:\n$self->{DIRECTORIES}{WORK}\n";
    print"Temporary files will be stored in\:\n$self->{DIRECTORIES}{TEMP}\n";
    print"Blast is located in\:\n$self->{DIRECTORIES}{BLAST}\n\n";
};

sub load_motif_set{
    my ($self, $set, $motifs) = @_;
    $self->{MOTIFS}{$set} = $motifs;
};

#Storage functions
sub store_genome{
    my $self = shift;
    my $file = shift;
    print "Storing genome in\:\n$file\n";
    my $hashref = \%{$self->{GENOME}};
    store $hashref, "$file";
    print "Done\n\n";
};

sub retrieve_genome{
    my $self = shift;
    my $file = shift;
    print "Loading genome from\:\n$file\n";
    my $holder = retrieve "$file";
    $self->{GENOME} = $holder;
    print "Done\n\n";
}

sub store_conserved{
    my $self = shift;
    my $file = shift;
    print "Storing conserved regions in\:\n$file\n";
    my $hashref = \%{$self->{CONSERVED}};
    store $hashref, "$file";
    print "Done\n\n";
};

sub retrieve_conserved{
    my $self = shift;
    my $file = shift;
    print "Loading conserved regions from\:\n$file\n";
    my $holder = retrieve "$file";
    $self->{CONSERVED} = $holder;
    print "Done\n\n";
}

sub store_data {
    my $self = shift;
    my $file = shift;
    print "Storing data in\:\n$file\n";
    store $self, "$file";
    print "Done\n\n";
};

sub retrieve_data {
    my $self = shift;
    my $file = shift;
    print "Loading annotation from\:\n$file\n";
    my $holder = retrieve "$file";
    $self->{SETS}   = $holder->{SETS};
    $self->{GENOME}   = $holder->{GENOME};
    $self->{CHROMOSOMES} = $holder->{CHROMOSOMES};
    $self->{CONSERVED} = $holder->{CONSERVED};
    $self->{NONCONSERVED} = $holder->{NONCONSERVED};
    $self->{DIRECTORIES} = $holder->{DIRECTORIES};
    $self->{INFO} = $holder->{INFO};
    my @test = keys %{$self->{GENOME}};
    print "Found ".($#test + 1)." records\n";
    print "Done\n\n";
};

sub load_overlaps {
    my $self = shift;
    my %holder = ();
    my %exon = ();
    my %intron = ();
    my @gene = keys %{$self->{GENOME}};
    print "Loading coding overlaps\n\n";
    for my $i (0 .. $#gene){
        for my $j (0 .. $#{$self->{GENOME}{$gene[$i]}{EXONSTART}}){
            my $record = { START => $self->{GENOME}{$gene[$i]}{EXONSTART}[$j],
                          END => $self->{GENOME}{$gene[$i]}{EXONEND}[$j]};
            push @{$holder{$self->{GENOME}{$gene[$i]}{CHROMOSOME}}}, $record;
        };
    };
    foreach my $Chromosome(keys %holder){
        @{$holder{$Chromosome}} = sort { $a->{START} <=> $b->{START} or $b->{END} <=> $a->{END}} @{$holder{$Chromosome}};
    };
    foreach my $Chromosome(keys %holder){
        push @{$exon{$Chromosome}}, $holder{$Chromosome}[0];
        for my $i ( 1 .. $#{$holder{$Chromosome}}){
            if ($holder{$Chromosome}[$i]->{END} < $exon{$Chromosome}[$#{$exon{$Chromosome}}]->{END}){next};
            if ($holder{$Chromosome}[$i]->{START} < $exon{$Chromosome}[$#{$exon{$Chromosome}}]->{END}){
                $exon{$Chromosome}[$#{$exon{$Chromosome}}]->{END} = $holder{$Chromosome}[$i]->{END}
            };
            if ($holder{$Chromosome}[$i]->{START} > $exon{$Chromosome}[$#{$exon{$Chromosome}}]->{END}){
                push @{$exon{$Chromosome}}, $holder{$Chromosome}[$i];
            };
        }
    };
    foreach my $Chromosome(keys %holder){
        my $record = {  START => '0',
                        END => $exon{$Chromosome}[0]->{START}};
        push @{$intron{$Chromosome}}, $record;
            
        for my $j (0 .. $#{$exon{$Chromosome}} - 1){
            my $element = {START => ($exon{$Chromosome}[$j]->{END} + 1 ),
                           END => ($exon{$Chromosome}[$j + 1]->{START} -1 )};
            push @{$intron{$Chromosome}}, $element;
        }
        my $last = {START => ($exon{$Chromosome}[$#{$exon{$Chromosome}}]->{END} + 1),
                    END => $self->end_chromosome($Chromosome)};
        push @{$intron{$Chromosome}}, $last;
    };
    $self->{EXONIC} = \%exon;
    $self->{INTRONIC} = \%intron;
};

#Gene extraction functions
sub list_upstream_primers {
    my ($self, $set, $maxlength, $primerlength, $searchlength, $directory) = @_;
    print "Listing upstream primers\n\n";
    my @gene = keys %{$self->{SETS}{$set}};
    @gene = sort { ($a =~ /(\d+)/)[0] <=> ($b =~ /(\d+)/)[0] } @gene;
    for my $i (0 .. $#gene){
        my ($leftbound, $rightbound) = undef;
        if ($self->{SETS}{$set}{$gene[$i]}{STRAND} =~ /\+/){
            $leftbound = $self->{SETS}{$set}{$gene[$i]}->{TXSTART} - $self->{INFO}{UPSTREAM};
            $rightbound = $self->{SETS}{$set}{$gene[$i]}->{TXSTART} + $self->{INFO}{DOWNSTREAM} - 1;
        }
        else{
            $leftbound = $self->{SETS}{$set}{$gene[$i]}->{TXEND} - $self->{INFO}{DOWNSTREAM} + 1;
            $rightbound = $self->{SETS}{$set}{$gene[$i]}->{TXEND} + $self->{INFO}{UPSTREAM};
        };
        my $chr = $self->{SETS}{$set}{$gene[$i]}->{CHROMOSOME};
        my @conserved = @{$self->get_conserved_blocks($chr, $leftbound, $rightbound)};
        my $count = 1;
        for my $j (0 .. $#conserved){
            my $width = $conserved[$j]->{END} - $conserved[$j]->{START};
            if ($width <= 0){next};
            my $frags = 1 + int($width/$maxlength);
            my $fraglength = int($width / $frags);
            for my $k (0 .. int($width/$maxlength)){
                my $start = $conserved[$j]->{START} + ($k * $fraglength);
                my $end = $conserved[$j]->{START} + (($k+1) * $fraglength);
                my $primer_start = $start - $searchlength;
                my $primer_end = $end + $searchlength;
                my $sequence = uppercase_sequence($self->get_sequence($chr, $primer_start, $primer_end));
                my $region_start = $searchlength + $primerlength;
                my $region_length = $fraglength - (2 * $primerlength) + 1;
                #print"$gene[$i]\_$count\nSequence\n$sequence\nStart\t$region_start\nLength\t$region_length\nSize\t".length($sequence)."\n\n";
                open (OUT, "> $directory/$gene[$i]\_$count\.pri");
                print OUT "PRIMER\_SEQUENCE\_ID\=$gene[$i]\_$count\n";
                print OUT "SEQUENCE\=$sequence\n";
                print OUT "TARGET\=$region_start\,$region_length\n";
                print OUT "PRIMER\_OPT\_SIZE\=$primerlength\n";
                print OUT "PRIMER\_MIN\_SIZE\=$primerlength\n";
                print OUT "PRIMER\_MAX\_SIZE\=$primerlength\n";
                print OUT "PRIMER\_OPT\_TM\=62\.0\n";
                print OUT "PRIMER\_MIN\_TM\=59\.0\n";
                print OUT "PRIMER\_MAX\_TM\=65\.0\n";
                print OUT "PRIMER\_PRODUCT\_SIZE\_RANGE\=$region_length\-".length($sequence)."\n";
                print OUT "PRIMER\_EXPLAIN_FLAG\=1\n";
                print OUT "PRIMER\_MISPRIMING_LIBRARY\=\n";
                print OUT "\=\n\n";
                close OUT;
                $count++;
            }
        }
    }
};

sub get_conserved_blocks {
    my $self = shift;
    my $chr = shift;
    my $startseq = shift;
    my $endseq = shift;
    my @holder = ();
    my ($startarray, $endarray) = $self->find_conserved_bounds($chr, $startseq, $endseq);
    for my $i ($startarray .. $endarray){
        my $maskstart = $self->{CONSERVED}{$chr}[$i]->{START};
        my $maskend = $self->{CONSERVED}{$chr}[$i]->{END};
        if ($self->{CONSERVED}{$chr}[$i]->{START} < $startseq){
            $maskstart = $startseq;
        };
        if ($self->{CONSERVED}{$chr}[$i]->{END} > $endseq){
            $maskend = $endseq; 
        };
        my $element = { START => $maskstart, END => $maskend };
        push @holder, $element;
    };
    return \@holder;
};

#Gene extraction functions
sub list_upstream {
    my ($self, $set, $file) = @_;    
    my @gene = keys %{$self->{SETS}{$set}};
    @gene = sort { ($a =~ /(\d+)/)[0] <=> ($b =~ /(\d+)/)[0] } @gene;
    print "Listing upstream regions\:\n$file\n\n";
    open(FILE, "> $file");
    for my $i (0 .. $#gene){
        my ($leftbound, $rightbound) = undef;
        if ($self->{SETS}{$set}{$gene[$i]}{STRAND} =~ /\+/){
            $leftbound = $self->{SETS}{$set}{$gene[$i]}->{TXSTART} - $self->{INFO}{UPSTREAM};
            $rightbound = $self->{SETS}{$set}{$gene[$i]}->{TXSTART} + $self->{INFO}{DOWNSTREAM} - 1;
        }
        else{
            $leftbound = $self->{SETS}{$set}{$gene[$i]}->{TXEND} - $self->{INFO}{DOWNSTREAM} + 1;
            $rightbound = $self->{SETS}{$set}{$gene[$i]}->{TXEND} + $self->{INFO}{UPSTREAM};
        };
        my $chr = $self->{SETS}{$set}{$gene[$i]}->{CHROMOSOME};
        my $sequence = $self->get_sequence($chr, $leftbound, $rightbound);
        $sequence = $self->mask_overlaps($chr, $leftbound, $rightbound, $gene[$i], $sequence);
        if ($self->{INFO}{CONSERVED}){
            $sequence = $self->mask_nonconserved($self->{SETS}{$set}{$gene[$i]}->{CHROMOSOME}, $leftbound, $rightbound, $sequence);
        };
        if ($self->{INFO}{MASKED}){$sequence = mask_repeats($sequence)};
        $self->{SETS}{$set}{$gene[$i]}->{SEQUENCES} = split_sequence($sequence, $leftbound);
        my $count = 1;
        if ($self->{SETS}{$set}{$gene[$i]}->{STRAND} =~ /\+/){
            for (my $j = 0; $j <= $#{$self->{SETS}{$set}{$gene[$i]}->{SEQUENCES}}; $j++){
                my $seq = $self->{SETS}{$set}{$gene[$i]}->{SEQUENCES}[$j]->{SEQUENCE};
                my $size = length($seq);
                if ($size < $self->{INFO}{MINIMUM}){next};
                $seq = uppercase_sequence($seq);
                print FILE "\>$gene[$i]\_$count\n$seq\n";
                $count++;
            };
        }
        else{
            for (my $j = $#{$self->{SETS}{$set}{$gene[$i]}->{SEQUENCES}}; $j >= 0; $j--){
                my $seq = $self->{SETS}{$set}{$gene[$i]}->{SEQUENCES}[$j]->{SEQUENCE};
                my $size = length($seq);
                if ($size < $self->{INFO}{MINIMUM}){next};
                $seq = reverse_complement($seq);
                $seq = uppercase_sequence($seq);
                print FILE "\>$gene[$i]\_$count\n$seq\n";
                $count++;
            };
        };
    };
    close FILE;
};

#Gene extraction functions
sub list_first_intron {
    my ($self, $set, $file) = @_;    
    my @gene = keys %{$self->{SETS}{$set}};
    @gene = sort { ($a =~ /(\d+)/)[0] <=> ($b =~ /(\d+)/)[0] } @gene;
    print "Listing first intron regions\:\n$file\n\n";
    open(FILE, "> $file");
    for my $i (0 .. $#gene){
        my ($leftbound, $rightbound) = undef;
        if (!exists $self->{SETS}{$set}{$gene[$i]}->{EXONSTART}[1]){next};
        if ($self->{SETS}{$set}{$gene[$i]}{STRAND} =~ /\+/){
            $leftbound = $self->{SETS}{$set}{$gene[$i]}->{EXONEND}[0] + 1;
            $rightbound = $self->{SETS}{$set}{$gene[$i]}->{EXONSTART}[1] - 1;
            if (($rightbound - $leftbound) > $self->{INFO}{UPSTREAM}){$rightbound = $leftbound + $self->{INFO}{UPSTREAM}};
        }
        else{
            $leftbound = $self->{SETS}{$set}{$gene[$i]}->{EXONEND}[$#{$self->{SETS}{$set}{$gene[$i]}->{EXONEND}} - 1] + 1;
            $rightbound = $self->{SETS}{$set}{$gene[$i]}->{EXONSTART}[$#{$self->{SETS}{$set}{$gene[$i]}->{EXONSTART}}] - 1;
            if (($rightbound - $leftbound) > $self->{INFO}{UPSTREAM}){$leftbound = $rightbound - $self->{INFO}{UPSTREAM}};
        };
        my $chr = $self->{SETS}{$set}{$gene[$i]}->{CHROMOSOME};
        my $sequence = $self->get_sequence($chr, $leftbound, $rightbound);
        $sequence = $self->mask_overlaps($chr, $leftbound, $rightbound, $gene[$i], $sequence);
        if ($self->{INFO}{CONSERVED}){
            $sequence = $self->mask_nonconserved($self->{SETS}{$set}{$gene[$i]}->{CHROMOSOME}, $leftbound, $rightbound, $sequence);
        };
        if ($self->{INFO}{MASKED}){$sequence = mask_repeats($sequence)};
        $self->{SETS}{$set}{$gene[$i]}->{SEQUENCES} = split_sequence($sequence, $leftbound);
        my $count = 1;
        if ($self->{SETS}{$set}{$gene[$i]}->{STRAND} =~ /\+/){
            for (my $j = 0; $j <= $#{$self->{SETS}{$set}{$gene[$i]}->{SEQUENCES}}; $j++){
                my $seq = $self->{SETS}{$set}{$gene[$i]}->{SEQUENCES}[$j]->{SEQUENCE};
                my $size = length($seq);
                if ($size < $self->{INFO}{MINIMUM}){next};
                $seq = uppercase_sequence($seq);
                print FILE "\>$gene[$i]\_$count\n$seq\n";
                $count++;
            };
        }
        else{
            for (my $j = $#{$self->{SETS}{$set}{$gene[$i]}->{SEQUENCES}}; $j >= 0; $j--){
                my $seq = $self->{SETS}{$set}{$gene[$i]}->{SEQUENCES}[$j]->{SEQUENCE};
                my $size = length($seq);
                if ($size < $self->{INFO}{MINIMUM}){next};
                $seq = reverse_complement($seq);
                $seq = uppercase_sequence($seq);
                print FILE "\>$gene[$i]\_$count\n$seq\n";
                $count++;
            };
        };
    };
    close FILE;
};

#Gene extraction functions
sub list_first_intron_primers {
    my ($self, $set, $maxlength, $primerlength, $searchlength, $directory) = @_;
    print "Listing first intron primers\n\n";
    my @gene = keys %{$self->{SETS}{$set}};
    @gene = sort { ($a =~ /(\d+)/)[0] <=> ($b =~ /(\d+)/)[0] } @gene;
    for my $i (0 .. $#gene){
        my ($leftbound, $rightbound) = undef;
        if (!exists $self->{SETS}{$set}{$gene[$i]}->{EXONSTART}[1]){next};
        if ($self->{SETS}{$set}{$gene[$i]}{STRAND} =~ /\+/){
            $leftbound = $self->{SETS}{$set}{$gene[$i]}->{EXONEND}[0] + 1;
            $rightbound = $self->{SETS}{$set}{$gene[$i]}->{EXONSTART}[1] - 1;
            if (($rightbound - $leftbound) > $self->{INFO}{UPSTREAM}){$rightbound = $leftbound + $self->{INFO}{UPSTREAM}};
        }
        else{
            $leftbound = $self->{SETS}{$set}{$gene[$i]}->{EXONEND}[$#{$self->{SETS}{$set}{$gene[$i]}->{EXONEND}} - 1] + 1;
            $rightbound = $self->{SETS}{$set}{$gene[$i]}->{EXONSTART}[$#{$self->{SETS}{$set}{$gene[$i]}->{EXONSTART}}] - 1;
            if (($rightbound - $leftbound) > $self->{INFO}{UPSTREAM}){$leftbound = $rightbound - $self->{INFO}{UPSTREAM}};
        };
        
        my $chr = $self->{SETS}{$set}{$gene[$i]}->{CHROMOSOME};
        my @conserved = @{$self->get_conserved_blocks($chr, $leftbound, $rightbound)};
        my $count = 1;
        for my $j (0 .. $#conserved){
            my $width = $conserved[$j]->{END} - $conserved[$j]->{START};
            if ($width <= 0){next};
            my $frags = 1 + int($width/$maxlength);
            my $fraglength = int($width / $frags);
            for my $k (0 .. int($width/$maxlength)){
                my $start = $conserved[$j]->{START} + ($k * $fraglength);
                my $end = $conserved[$j]->{START} + (($k+1) * $fraglength);
                my $primer_start = $start - $searchlength;
                my $primer_end = $end + $searchlength;
                my $sequence = uppercase_sequence($self->get_sequence($chr, $primer_start, $primer_end));
                my $region_start = $searchlength + $primerlength;
                my $region_length = $fraglength - (2 * $primerlength) + 1;
                #print"$gene[$i]\_$count\nSequence\n$sequence\nStart\t$region_start\nLength\t$region_length\nSize\t".length($sequence)."\n\n";
                open (OUT, "> $directory/$gene[$i]\_$count\.pri");
                print OUT "PRIMER\_SEQUENCE\_ID\=$gene[$i]\_$count\n";
                print OUT "SEQUENCE\=$sequence\n";
                print OUT "TARGET\=$region_start\,$region_length\n";
                print OUT "PRIMER\_OPT\_SIZE\=$primerlength\n";
                print OUT "PRIMER\_MIN\_SIZE\=$primerlength\n";
                print OUT "PRIMER\_MAX\_SIZE\=$primerlength\n";
                print OUT "PRIMER\_OPT\_TM\=62\.0\n";
                print OUT "PRIMER\_MIN\_TM\=59\.0\n";
                print OUT "PRIMER\_MAX\_TM\=65\.0\n";
                print OUT "PRIMER\_PRODUCT\_SIZE\_RANGE\=$region_length\-".length($sequence)."\n";
                print OUT "PRIMER\_EXPLAIN_FLAG\=1\n";
                print OUT "PRIMER\_MISPRIMING_LIBRARY\=\n";
                print OUT "\=\n\n";
                close OUT;
                $count++;
            }
        }
    }
};

#Gene extraction functions
sub list_downstream {
    my ($self, $set, $file) = @_;    
    my @gene = keys %{$self->{SETS}{$set}};
    @gene = sort { ($a =~ /(\d+)/)[0] <=> ($b =~ /(\d+)/)[0] } @gene;
    print "Listing downstream regions\:\n$file\n\n";
    open(FILE, "> $file");
    for my $i (0 .. $#gene){
        my ($leftbound, $rightbound) = undef;
        if ($self->{SETS}{$set}{$gene[$i]}{STRAND} =~ /\+/){
            $leftbound = $self->{SETS}{$set}{$gene[$i]}->{TXEND} - $self->{INFO}{DOWNSTREAM} + 1;
            $rightbound = $self->{SETS}{$set}{$gene[$i]}->{TXEND} + $self->{INFO}{UPSTREAM};
        }
        else{
            $leftbound = $self->{SETS}{$set}{$gene[$i]}->{TXSTART} - $self->{INFO}{UPSTREAM};
            $rightbound = $self->{SETS}{$set}{$gene[$i]}->{TXSTART} + $self->{INFO}{DOWNSTREAM} - 1;
        };
        my $chr = $self->{SETS}{$set}{$gene[$i]}->{CHROMOSOME};
        my $sequence = $self->get_sequence($chr, $leftbound, $rightbound);
        $sequence = $self->mask_overlaps($chr, $leftbound, $rightbound, $gene[$i], $sequence);
        if ($self->{INFO}{CONSERVED}){
            $sequence = $self->mask_nonconserved($self->{SETS}{$set}{$gene[$i]}->{CHROMOSOME}, $leftbound, $rightbound, $sequence);
        };
        if ($self->{INFO}{MASKED}){$sequence = mask_repeats($sequence)};
        $self->{SETS}{$set}{$gene[$i]}->{SEQUENCES} = split_sequence($sequence, $leftbound);
        my $count = 1;
        if ($self->{SETS}{$set}{$gene[$i]}->{STRAND} =~ /\+/){
            for (my $j = 0; $j <= $#{$self->{SETS}{$set}{$gene[$i]}->{SEQUENCES}}; $j++){
                my $seq = $self->{SETS}{$set}{$gene[$i]}->{SEQUENCES}[$j]->{SEQUENCE};
                my $size = length($seq);
                if ($size < $self->{INFO}{MINIMUM}){next};
                $seq = uppercase_sequence($seq);
                print FILE "\>$gene[$i]\_$count\n$seq\n";
                $count++;
            };
        }
        else{
            for (my $j = $#{$self->{SETS}{$set}{$gene[$i]}->{SEQUENCES}}; $j >= 0; $j--){
                my $seq = $self->{SETS}{$set}{$gene[$i]}->{SEQUENCES}[$j]->{SEQUENCE};
                my $size = length($seq);
                if ($size < $self->{INFO}{MINIMUM}){next};
                $seq = reverse_complement($seq);
                $seq = uppercase_sequence($seq);
                print FILE "\>$gene[$i]\_$count\n$seq\n";
                $count++;
            };
        };
    };
    close FILE;
};

#Gene extraction functions
sub list_downstream_primers {
    my ($self, $set, $maxlength, $primerlength, $searchlength, $directory) = @_;    
    my @gene = keys %{$self->{SETS}{$set}};
    print "Listing downstream primers\n\n";
    @gene = sort { ($a =~ /(\d+)/)[0] <=> ($b =~ /(\d+)/)[0] } @gene;
    for my $i (0 .. $#gene){
        my ($leftbound, $rightbound) = undef;
        if ($self->{SETS}{$set}{$gene[$i]}{STRAND} =~ /\+/){
            $leftbound = $self->{SETS}{$set}{$gene[$i]}->{TXEND} - $self->{INFO}{DOWNSTREAM} + 1;
            $rightbound = $self->{SETS}{$set}{$gene[$i]}->{TXEND} + $self->{INFO}{UPSTREAM};
        }
        else{
            $leftbound = $self->{SETS}{$set}{$gene[$i]}->{TXSTART} - $self->{INFO}{UPSTREAM};
            $rightbound = $self->{SETS}{$set}{$gene[$i]}->{TXSTART} + $self->{INFO}{DOWNSTREAM} - 1;
        };
        my $chr = $self->{SETS}{$set}{$gene[$i]}->{CHROMOSOME};
        my @conserved = @{$self->get_conserved_blocks($chr, $leftbound, $rightbound)};
        my $count = 1;
        for my $j (0 .. $#conserved){
            my $width = $conserved[$j]->{END} - $conserved[$j]->{START};
            if ($width <= 0){next};
            my $frags = 1 + int($width/$maxlength);
            my $fraglength = int($width / $frags);
            for my $k (0 .. int($width/$maxlength)){
                my $start = $conserved[$j]->{START} + ($k * $fraglength);
                my $end = $conserved[$j]->{START} + (($k+1) * $fraglength);
                my $primer_start = $start - $searchlength;
                my $primer_end = $end + $searchlength;
                my $sequence = uppercase_sequence($self->get_sequence($chr, $primer_start, $primer_end));
                my $region_start = $searchlength + $primerlength;
                my $region_length = $fraglength - (2 * $primerlength) + 1;
                #print"$gene[$i]\_$count\nSequence\n$sequence\nStart\t$region_start\nLength\t$region_length\nSize\t".length($sequence)."\n\n";
                open (OUT, "> $directory/$gene[$i]\_$count\.pri");
                print OUT "PRIMER\_SEQUENCE\_ID\=$gene[$i]\_$count\n";
                print OUT "SEQUENCE\=$sequence\n";
                print OUT "TARGET\=$region_start\,$region_length\n";
                print OUT "PRIMER\_OPT\_SIZE\=$primerlength\n";
                print OUT "PRIMER\_MIN\_SIZE\=$primerlength\n";
                print OUT "PRIMER\_MAX\_SIZE\=$primerlength\n";
                print OUT "PRIMER\_OPT\_TM\=62\.0\n";
                print OUT "PRIMER\_MIN\_TM\=59\.0\n";
                print OUT "PRIMER\_MAX\_TM\=65\.0\n";
                print OUT "PRIMER\_PRODUCT\_SIZE\_RANGE\=$region_length\-".length($sequence)."\n";
                print OUT "PRIMER\_EXPLAIN_FLAG\=1\n";
                print OUT "PRIMER\_MISPRIMING_LIBRARY\=\n";
                print OUT "\=\n\n";
                close OUT;
                $count++;
            }
        }
    }
};

#Gene extraction functions
sub list_genome_upstream {
    my ($self, $file) = @_;    
    my @gene = keys %{$self->{GENOME}};
    @gene = sort { ($a =~ /(\d+)/)[0] <=> ($b =~ /(\d+)/)[0] } @gene;
    print "Listing genome sequences\:\n$file\n\n";
    open(FILE, "> $file");
    for my $i (0 .. $#gene){
        my ($leftbound, $rightbound) = undef;
        if ($self->{GENOME}{$gene[$i]}{STRAND} =~ /\+/){
            $leftbound = $self->{GENOME}{$gene[$i]}->{TXSTART} - $self->{INFO}{UPSTREAM};
            $rightbound = $self->{GENOME}{$gene[$i]}->{TXSTART} + $self->{INFO}{DOWNSTREAM} - 1;
        }
        else{
            $leftbound = $self->{GENOME}{$gene[$i]}->{TXEND} - $self->{INFO}{DOWNSTREAM} + 1;
            $rightbound = $self->{GENOME}{$gene[$i]}->{TXEND} + $self->{INFO}{UPSTREAM};
        };
        my $chr = $self->{GENOME}{$gene[$i]}->{CHROMOSOME};
        my $sequence = $self->get_sequence($chr, $leftbound, $rightbound);
        $sequence = $self->mask_overlaps($chr, $leftbound, $rightbound, $gene[$i], $sequence);
        if ($self->{INFO}{CONSERVED}){
            $sequence = $self->mask_nonconserved($self->{GENOME}{$gene[$i]}->{CHROMOSOME}, $leftbound, $rightbound, $sequence);
        };
        if ($self->{INFO}{MASKED}){$sequence = mask_repeats($sequence)};
        $self->{GENOME}{$gene[$i]}->{SEQUENCES} = split_sequence($sequence, $leftbound);
        my $count = 1;
        if ($self->{GENOME}{$gene[$i]}->{STRAND} =~ /\+/){
            for (my $j = 0; $j <= $#{$self->{GENOME}{$gene[$i]}->{SEQUENCES}}; $j++){
                my $seq = $self->{GENOME}{$gene[$i]}->{SEQUENCES}[$j]->{SEQUENCE};
                my $size = length($seq);
                if ($size < $self->{INFO}{MINIMUM}){next};
                $seq = uppercase_sequence($seq);
                print FILE "\>$gene[$i]\_$count\n$seq\n";
                $count++;
            };
        }
        else{
            for (my $j = $#{$self->{GENOME}{$gene[$i]}->{SEQUENCES}}; $j >= 0; $j--){
                my $seq = $self->{GENOME}{$gene[$i]}->{SEQUENCES}[$j]->{SEQUENCE};
                my $size = length($seq);
                if ($size < $self->{INFO}{MINIMUM}){next};
                $seq = reverse_complement($seq);
                $seq = uppercase_sequence($seq);
                print FILE "\>$gene[$i]\_$count\n$seq\n";
                $count++;
            };
        };
    };
    close FILE;
};

#Gene extraction functions
sub list_genome_first_intron {
    my ($self, $file) = @_;    
    my @gene = keys %{$self->{GENOME}};
    @gene = sort { ($a =~ /(\d+)/)[0] <=> ($b =~ /(\d+)/)[0] } @gene;
    print "Listing genome intron sequences\:\n$file\n\n";
    open(FILE, "> $file");
    for my $i (0 .. $#gene){
        my ($leftbound, $rightbound) = undef;
        if (!exists $self->{GENOME}{$gene[$i]}->{EXONSTART}[1]){next};
        if ($self->{GENOME}{$gene[$i]}{STRAND} =~ /\+/){
            $leftbound = $self->{GENOME}{$gene[$i]}->{EXONEND}[0] + 1;
            $rightbound = $self->{GENOME}{$gene[$i]}->{EXONSTART}[1] - 1;
            if (($rightbound - $leftbound) > $self->{INFO}{UPSTREAM}){$rightbound = $leftbound + $self->{INFO}{UPSTREAM}};
        }
        else{
            $leftbound = $self->{GENOME}{$gene[$i]}->{EXONEND}[$#{$self->{GENOME}{$gene[$i]}->{EXONEND}} - 1] + 1;
            $rightbound = $self->{GENOME}{$gene[$i]}->{EXONSTART}[$#{$self->{GENOME}{$gene[$i]}->{EXONSTART}}] - 1;
            if (($rightbound - $leftbound) > $self->{INFO}{UPSTREAM}){$leftbound = $rightbound - $self->{INFO}{UPSTREAM}};
        };
        my $chr = $self->{GENOME}{$gene[$i]}->{CHROMOSOME};
        my $sequence = $self->get_sequence($chr, $leftbound, $rightbound);
        $sequence = $self->mask_overlaps($chr, $leftbound, $rightbound, $gene[$i], $sequence);
        if ($self->{INFO}{CONSERVED}){
            $sequence = $self->mask_nonconserved($self->{GENOME}{$gene[$i]}->{CHROMOSOME}, $leftbound, $rightbound, $sequence);
        };
        if ($self->{INFO}{MASKED}){$sequence = mask_repeats($sequence)};
        $self->{GENOME}{$gene[$i]}->{SEQUENCES} = split_sequence($sequence, $leftbound);
        
        my $count = 1;
        if ($self->{GENOME}{$gene[$i]}->{STRAND} =~ /\+/){
            for (my $j = 0; $j <= $#{$self->{GENOME}{$gene[$i]}->{SEQUENCES}}; $j++){
                my $seq = $self->{GENOME}{$gene[$i]}->{SEQUENCES}[$j]->{SEQUENCE};
                my $size = length($seq);
                if ($size < $self->{INFO}{MINIMUM}){next};
                $seq = uppercase_sequence($seq);
                print FILE "\>$gene[$i]\_$count\n$seq\n";
                $count++;
            };
        }
        else{
            for (my $j = $#{$self->{GENOME}{$gene[$i]}->{SEQUENCES}}; $j >= 0; $j--){
                my $seq = $self->{GENOME}{$gene[$i]}->{SEQUENCES}[$j]->{SEQUENCE};
                my $size = length($seq);
                if ($size < $self->{INFO}{MINIMUM}){next};
                $seq = reverse_complement($seq);
                $seq = uppercase_sequence($seq);
                print FILE "\>$gene[$i]\_$count\n$seq\n";
                $count++;
            };
        };
    };
    close FILE;
};

#Gene extraction functions
sub list_genome_downstream {
    my ($self, $file) = @_;    
    my @gene = keys %{$self->{GENOME}};
    @gene = sort { ($a =~ /(\d+)/)[0] <=> ($b =~ /(\d+)/)[0] } @gene;
    print "Listing genome sequences\:\n$file\n\n";
    open(FILE, "> $file");
    for my $i (0 .. $#gene){
        my ($leftbound, $rightbound) = undef;
        if (!defined $self->{GENOME}{$gene[$i]}{STRAND}){next};
        if ($self->{GENOME}{$gene[$i]}{STRAND} =~ /\+/){
            $leftbound = $self->{GENOME}{$gene[$i]}->{TXEND} - $self->{INFO}{DOWNSTREAM} + 1;
            $rightbound = $self->{GENOME}{$gene[$i]}->{TXEND} + $self->{INFO}{UPSTREAM};
        }
        else{
            $leftbound = $self->{GENOME}{$gene[$i]}->{TXSTART} - $self->{INFO}{UPSTREAM};
            $rightbound = $self->{GENOME}{$gene[$i]}->{TXSTART} + $self->{INFO}{DOWNSTREAM} - 1;
        };
        my $chr = $self->{GENOME}{$gene[$i]}->{CHROMOSOME};
        my $sequence = $self->get_sequence($chr, $leftbound, $rightbound);
        $sequence = $self->mask_overlaps($chr, $leftbound, $rightbound, $gene[$i], $sequence);
        if ($self->{INFO}{CONSERVED}){
            $sequence = $self->mask_nonconserved($self->{GENOME}{$gene[$i]}->{CHROMOSOME}, $leftbound, $rightbound, $sequence);
        };
        if ($self->{INFO}{MASKED}){$sequence = mask_repeats($sequence)};
        $self->{GENOME}{$gene[$i]}->{SEQUENCES} = split_sequence($sequence, $leftbound);
        my $count = 1;
        if ($self->{GENOME}{$gene[$i]}->{STRAND} =~ /\+/){
            for (my $j = 0; $j <= $#{$self->{GENOME}{$gene[$i]}->{SEQUENCES}}; $j++){
                my $seq = $self->{GENOME}{$gene[$i]}->{SEQUENCES}[$j]->{SEQUENCE};
                my $size = length($seq);
                if ($size < $self->{INFO}{MINIMUM}){next};
                $seq = uppercase_sequence($seq);
                print FILE "\>$gene[$i]\_$count\n$seq\n";
                $count++;
            };
        }
        else{
            for (my $j = $#{$self->{GENOME}{$gene[$i]}->{SEQUENCES}}; $j >= 0; $j--){
                my $seq = $self->{GENOME}{$gene[$i]}->{SEQUENCES}[$j]->{SEQUENCE};
                my $size = length($seq);
                if ($size < $self->{INFO}{MINIMUM}){next};
                $seq = reverse_complement($seq);
                $seq = uppercase_sequence($seq);
                print FILE "\>$gene[$i]\_$count\n$seq\n";
                $count++;
            };
        };
    };
    close FILE;
};

#Gene extraction functions
sub list_scanace_upstream {
    my ($self, $file) = @_;    
    my @gene = keys %{$self->{GENOME}};
    @gene = sort { ($a =~ /(\d+)/)[0] <=> ($b =~ /(\d+)/)[0] } @gene;
    print "Listing scanace sequences\:\n$file\n\n";
    open(FILE, "> $file");
    for my $i (0 .. $#gene){
        my ($leftbound, $rightbound) = undef;
        if ($self->{GENOME}{$gene[$i]}{STRAND} =~ /\+/){
            $leftbound = $self->{GENOME}{$gene[$i]}->{TXSTART} - $self->{INFO}{UPSTREAM};
            $rightbound = $self->{GENOME}{$gene[$i]}->{TXSTART} + $self->{INFO}{DOWNSTREAM} - 1;
        }
        else{
            $leftbound = $self->{GENOME}{$gene[$i]}->{TXEND} - $self->{INFO}{DOWNSTREAM} + 1;
            $rightbound = $self->{GENOME}{$gene[$i]}->{TXEND} + $self->{INFO}{UPSTREAM};
        };
        my $chr = $self->{GENOME}{$gene[$i]}->{CHROMOSOME};
        my $sequence = $self->get_sequence($chr, $leftbound, $rightbound);
        $sequence = $self->mask_overlaps($chr, $leftbound, $rightbound, $gene[$i], $sequence);
        if ($self->{INFO}{CONSERVED}){
            $sequence = $self->mask_nonconserved($self->{GENOME}{$gene[$i]}->{CHROMOSOME}, $leftbound, $rightbound, $sequence);
        };
        if ($self->{INFO}{MASKED}){$sequence = mask_repeats($sequence)};
        if ($self->{GENOME}{$gene[$i]}->{STRAND} =~ /\-/){$sequence = reverse_complement($sequence)};
        $sequence = uppercase_sequence($sequence);
        print FILE "\>$gene[$i]\n$sequence\n";
    };
    close FILE;
};

#Gene extraction functions
sub list_scanace_first_intron {
    my ($self, $file) = @_;    
    my @gene = keys %{$self->{GENOME}};
    @gene = sort { ($a =~ /(\d+)/)[0] <=> ($b =~ /(\d+)/)[0] } @gene;
    print "Listing scanace sequences\:\n$file\n\n";
    open(FILE, "> $file");
    for my $i (0 .. $#gene){
        my ($leftbound, $rightbound) = undef;
        if (!exists $self->{GENOME}{$gene[$i]}->{EXONSTART}[1]){next};
        if ($self->{GENOME}{$gene[$i]}{STRAND} =~ /\+/){
            $leftbound = $self->{GENOME}{$gene[$i]}->{EXONEND}[0] + 1;
            $rightbound = $self->{GENOME}{$gene[$i]}->{EXONSTART}[1] - 1;
            if (($rightbound - $leftbound) > $self->{INFO}{UPSTREAM}){$rightbound = $leftbound + $self->{INFO}{UPSTREAM}};
        }
        else{
            $leftbound = $self->{GENOME}{$gene[$i]}->{EXONEND}[$#{$self->{GENOME}{$gene[$i]}->{EXONEND}} - 1] + 1;
            $rightbound = $self->{GENOME}{$gene[$i]}->{EXONSTART}[$#{$self->{GENOME}{$gene[$i]}->{EXONSTART}}] - 1;
            if (($rightbound - $leftbound) > $self->{INFO}{UPSTREAM}){$leftbound = $rightbound - $self->{INFO}{UPSTREAM}};
        };
        my $chr = $self->{GENOME}{$gene[$i]}->{CHROMOSOME};
        my $sequence = $self->get_sequence($chr, $leftbound, $rightbound);
        $sequence = $self->mask_overlaps($chr, $leftbound, $rightbound, $gene[$i], $sequence);
        if ($self->{INFO}{CONSERVED}){
            $sequence = $self->mask_nonconserved($self->{GENOME}{$gene[$i]}->{CHROMOSOME}, $leftbound, $rightbound, $sequence);
        };
        if ($self->{INFO}{MASKED}){$sequence = mask_repeats($sequence)};
        if ($self->{GENOME}{$gene[$i]}->{STRAND} =~ /\-/){$sequence = reverse_complement($sequence)};
        $sequence = uppercase_sequence($sequence);
        print FILE "\>$gene[$i]\n$sequence\n";
    };
    close FILE;
};

#Gene extraction functions
sub list_scanace_downstream {
    my ($self, $file) = @_;    
    my @gene = keys %{$self->{GENOME}};
    @gene = sort { ($a =~ /(\d+)/)[0] <=> ($b =~ /(\d+)/)[0] } @gene;
    print "Listing scanace sequences\:\n$file\n\n";
    open(FILE, "> $file");
    for my $i (0 .. $#gene){
        my ($leftbound, $rightbound) = undef;
        if (!defined $self->{GENOME}{$gene[$i]}{STRAND}){next};
        if ($self->{GENOME}{$gene[$i]}{STRAND} =~ /\+/){
            $leftbound = $self->{GENOME}{$gene[$i]}->{TXEND} - $self->{INFO}{DOWNSTREAM} + 1;
            $rightbound = $self->{GENOME}{$gene[$i]}->{TXEND} + $self->{INFO}{UPSTREAM};
        }
        else{
            $leftbound = $self->{GENOME}{$gene[$i]}->{TXSTART} - $self->{INFO}{UPSTREAM};
            $rightbound = $self->{GENOME}{$gene[$i]}->{TXSTART} + $self->{INFO}{DOWNSTREAM} - 1;
        };
        my $chr = $self->{GENOME}{$gene[$i]}->{CHROMOSOME};
        my $sequence = $self->get_sequence($chr, $leftbound, $rightbound);
        $sequence = $self->mask_overlaps($chr, $leftbound, $rightbound, $gene[$i], $sequence);
        if ($self->{INFO}{CONSERVED}){
            $sequence = $self->mask_nonconserved($self->{GENOME}{$gene[$i]}->{CHROMOSOME}, $leftbound, $rightbound, $sequence);
        };
        if ($self->{INFO}{MASKED}){$sequence = mask_repeats($sequence)};
        if ($self->{GENOME}{$gene[$i]}->{STRAND} =~ /\-/){$sequence = reverse_complement($sequence)};
        $sequence = uppercase_sequence($sequence);
        print FILE "\>$gene[$i]\n$sequence\n";
    };
    close FILE;
};

sub mask_overlaps {
    my $self = shift;
    my $chr = shift;
    my $startseq = shift;
    my $endseq = shift;
    my $gene = shift;
    my $sequence = shift;
    if (defined $self->{GENOME}{$gene}->{OVERLAPS}){
    my @overlaps = keys %{$self->{GENOME}{$gene}->{OVERLAPS}};
    for my $i (0 .. $#overlaps){
        for my $j (0 .. $#{$self->{GENOME}{$gene}->{OVERLAPS}{$overlaps[$i]}}){
            my $maskstart = $self->{GENOME}{$gene}->{OVERLAPS}{$overlaps[$i]}[$j]->{START};
            my $maskend = $self->{GENOME}{$gene}->{OVERLAPS}{$overlaps[$i]}[$j]->{END};
            if ($maskend < $startseq){next};
            if ($maskstart > $endseq){next};
            if ($self->{GENOME}{$gene}->{OVERLAPS}{$overlaps[$i]}[$j]->{START} < $startseq){
                $maskstart = $startseq;
            };
            if ($self->{GENOME}{$gene}->{OVERLAPS}{$overlaps[$i]}[$j]->{END} > $endseq){
                $maskend = $endseq; 
            };
            my $beginmask =  $maskstart - $startseq;
            my $length = $maskend - $maskstart + 1;
            if ($length < 0){next};
            my $mask = 'N' x $length;
            substr($sequence, $beginmask, $length, $mask);
        };
    };
    };
    return $sequence;
};

#sequence functions
sub get_sequence {
    my $self = shift;
    my $chr = shift;
    my $startseq = shift;
    my $endseq = shift;
    my $sequence = '';
    my $chromfile = "$self->{DIRECTORIES}{GENOME}\/chr$chr\.fa";
    $/ = "\n";
    open(SEQ, "< $chromfile");
    $. = 0; #finds the length of the first line
    my $first_line = <SEQ>;
    my $length_first_line = length($first_line);
    $. = 1; #finds the length of the second line and the record separator
    my $second_line = <SEQ>;
    chomp $second_line;
    my $length_second_line = length ($second_line);
    if ($length_second_line < 1){
        print "Problem with this local $chr\t$startseq\t$endseq\n";
        return $sequence};
    my $start = $startseq + int($startseq/$length_second_line) + $length_first_line - 1;
    my $end = $endseq + int($endseq/$length_second_line) + $length_first_line - 1;
    my $length = $end - $start + 1;
    if ($length <= 0){return($sequence)};
    seek(SEQ, $start, 0);
    read(SEQ, $sequence, $length);
    close SEQ;
    $sequence =~ s/\n//g;
    return $sequence;
};

sub mask_repeats {
    my $sequence = shift;
    $sequence =~ tr/agtcn/NNNNN/;
    return $sequence;
};

sub reverse_complement {
    my $seq = shift;
    $seq =~ tr/ACTGactg/TGACtgac/;
    $seq = reverse $seq;
    return $seq;
};

sub uppercase_sequence {
    my $sequence = shift;
    $sequence =~ tr/agtcn/AGTCN/;
    return $sequence;
};

sub split_nonrepeat{
    my $sequence = shift;
    my $start = shift;
    $sequence =~ s/([AGTC])(N)/$1\t$2/g;
    $sequence =~ s/(N)([AGTC])/$1\t$2/g;
    my @fragment = split(/\t/, $sequence);
    my @nonrepeats = ();
    for my $z (0 .. $#fragment){
        if ($fragment[$z] =~ /N+/) {next};
        if (length($fragment[$z]) < 1){next};
        my $frag_start = $start;
        if ($frag_start < 1){next};
        my $frag_end = $start -1 ;
        if ($z > 0){for my $a (0 .. $z-1){$frag_start += length($fragment[$a])}};
        for my $a (0 .. $z){$frag_end += length($fragment[$a])};
        my $element = {START => $frag_start, END => $frag_end};
        push @nonrepeats, $element;
    };
    return (\@nonrepeats);
}

sub split_sequence{
    my $sequence = shift;
    my $start = shift;
    $sequence =~ s/([AGTC])(N)/$1\t$2/g;
    $sequence =~ s/(N)([AGTC])/$1\t$2/g;
    my @fragment = split(/\t/, $sequence);
    my @seqs = ();
    for my $z (0 .. $#fragment){
        if ($fragment[$z] =~ /N+/) {next};
        if (length($fragment[$z]) < 1){next};
        my $frag_start = $start;
        if ($frag_start < 1){next};
        my $frag_end = $start - 1;
        if ($z > 0){for my $a (0 .. $z-1){$frag_start += length($fragment[$a])}};
        for my $a (0 .. $z){$frag_end += length($fragment[$a])};
        my $element = {START => $frag_start, END => $frag_end, SEQUENCE => $fragment[$z]};
        push @seqs, $element;
    };
    return (\@seqs);
}

sub find_nonrepeat{
    my ($self, $set) = @_;
    my @gene = keys %{$self->{SETS}{$set}};
    @gene = sort { ($a =~ /(\d+)/)[0] <=> ($b =~ /(\d+)/)[0] }  @gene;
    print"Finding nonrepeat regions\n";
    my $count = 0;
    for my $i (0 .. $#gene){
        my ($leftbound, $rightbound) = undef;
        if ($self->{SETS}{$set}{$gene[$i]}{STRAND} =~ /\+/){
            $leftbound = $self->{SETS}{$set}{$gene[$i]}->{TXSTART} - $self->{INFO}{UPSTREAM};
            $rightbound = $self->{SETS}{$set}{$gene[$i]}->{TXEND} + $self->{INFO}{UPSTREAM} - 1;
        }
        else{
            $leftbound = $self->{SETS}{$set}{$gene[$i]}->{TXSTART} - $self->{INFO}{UPSTREAM} + 1;
            $rightbound = $self->{SETS}{$set}{$gene[$i]}->{TXEND} + $self->{INFO}{UPSTREAM};
        };
        my $chr = $self->{SETS}{$set}{$gene[$i]}->{CHROMOSOME};
        my $sequence = $self->get_sequence($chr, $leftbound, $rightbound);
        $sequence = mask_repeats($sequence);
        $self->{SETS}{$set}{$gene[$i]}->{NONREPEAT} = split_nonrepeat($sequence, $leftbound);
        if ($count > 1000){print'.';$count = 0};
        $count++;
    };
    print"\n\n";
};

sub mask_nonconserved {
    my $self = shift;
    my $chr = shift;
    my $startseq = shift;
    my $endseq = shift;
    my $sequence = shift;
    my ($startarray, $endarray) = $self->find_nonconserved_bounds($chr, $startseq, $endseq);
    for my $i ($startarray .. $endarray){
        my $maskstart = $self->{NONCONSERVED}{$chr}[$i]->{START};
        my $maskend = $self->{NONCONSERVED}{$chr}[$i]->{END};
        if ($self->{NONCONSERVED}{$chr}[$i]->{START} < $startseq){
            $maskstart = $startseq;
        };
        if ($self->{NONCONSERVED}{$chr}[$i]->{END} > $endseq){
            $maskend = $endseq; 
        };
        my $beginmask =  $maskstart - $startseq;
        my $length = $maskend - $maskstart + 1;
        if ($length < 0){next};
        my $mask = 'N' x $length;
        substr($sequence, $beginmask, $length, $mask); 
    };
    return $sequence;
};

sub find_nonconserved_bounds{
    my $self = shift;
    my $chr = shift;
    my $leftbound = shift;
    my $rightbound = shift;
    my $leftarraybound = 0;
    my $rightarraybound = $#{$self->{NONCONSERVED}{$chr}};
    my ($leftcheck, $rightcheck) = undef;
    my $middle = int(($rightarraybound + $leftarraybound)/2);
    my $check = 1;
    my $loopcounter = 0;
    while ($check){
        $leftcheck = $leftarraybound;
        $rightcheck = $rightarraybound;
        if ($self->{NONCONSERVED}{$chr}[$middle]->{START} < $leftbound){
            $leftarraybound = $middle;
        };
        if ($self->{NONCONSERVED}{$chr}[$middle]->{END} > $rightbound){
            $rightarraybound = $middle;
        };
        $middle = int(($rightarraybound + $leftarraybound)/2);
        if ($leftarraybound == $leftcheck and $rightarraybound == $rightcheck){$check = 0};
        if ($loopcounter > 100){$check = 0};
        $loopcounter++;
    };
    return ($leftarraybound, $rightarraybound);
};

sub find_conserved_bounds{
    my $self = shift;
    my $chr = shift;
    my $leftbound = shift;
    my $rightbound = shift;
    my $leftarraybound = 0;
    my $rightarraybound = $#{$self->{CONSERVED}{$chr}};
    my ($leftcheck, $rightcheck) = undef;
    my $middle = int(($rightarraybound + $leftarraybound)/2);
    my $check = 1;
    my $loopcounter = 0;
    while ($check){
        $leftcheck = $leftarraybound;
        $rightcheck = $rightarraybound;
        if ($self->{CONSERVED}{$chr}[$middle]->{START} < $leftbound){
            $leftarraybound = $middle;
        };
        if ($self->{CONSERVED}{$chr}[$middle]->{END} > $rightbound){
            $rightarraybound = $middle;
        };
        $middle = int(($rightarraybound + $leftarraybound)/2);
        if ($leftarraybound == $leftcheck and $rightarraybound == $rightcheck){$check = 0};
        if ($loopcounter > 100){$check = 0};
        $loopcounter++;
    };
    return ($leftarraybound, $rightarraybound);
};

sub find_factor_bounds{
    my $self = shift;
    my $chr = shift;
    my $leftbound = shift;
    my $rightbound = shift;
    my $leftarraybound = 0;
    my $rightarraybound = $#{$self->{FACTORS}{$chr}};
    my ($leftcheck, $rightcheck) = undef;
    my $middle = int(($rightarraybound + $leftarraybound)/2);
    my $check = 1;
    my $loopcounter = 0;
    while ($check){
        $leftcheck = $leftarraybound;
        $rightcheck = $rightarraybound;
        if ($self->{FACTORS}{$chr}[$middle]->{START} < $leftbound){
            $leftarraybound = $middle;
        };
        if ($self->{FACTORS}{$chr}[$middle]->{END} > $rightbound){
            $rightarraybound = $middle;
        };
        $middle = int(($rightarraybound + $leftarraybound)/2);
        if ($leftarraybound == $leftcheck and $rightarraybound == $rightcheck){$check = 0};
        if ($loopcounter > 100){$check = 0};
        $loopcounter++;
    };
    return ($leftarraybound, $rightarraybound);
};

sub find_nonconserved{
    my $self = shift;
    my @gene = keys %{$self->{GENOME}};
    print"Finding nonconserved regions\n";
    my $count = 0;
    for my $i (0 .. $#gene){
        my ($leftbound, $rightbound) = undef;
        if ($self->{GENOME}{$gene[$i]}{STRAND} =~ /\+/){
            $leftbound = $self->{GENOME}{$gene[$i]}->{TXSTART} - $self->{INFO}{UPSTREAM};
            $rightbound = $self->{GENOME}{$gene[$i]}->{TXEND} + $self->{INFO}{UPSTREAM} - 1;
        }
        else{
            $leftbound = $self->{GENOME}{$gene[$i]}->{TXSTART} - $self->{INFO}{UPSTREAM} + 1;
            $rightbound = $self->{GENOME}{$gene[$i]}->{TXEND} + $self->{INFO}{UPSTREAM};
        };
        my $chr = $self->{GENOME}{$gene[$i]}->{CHROMOSOME};
        my ($up, $down) = $self->find_nonconserved_bounds($chr, $leftbound, $rightbound);
        for my $j ($up .. $down){
            my $maskstart = $self->{NONCONSERVED}{$chr}[$j]->{START};
            my $maskend = $self->{NONCONSERVED}{$chr}[$j]->{END};
            if ($self->{NONCONSERVED}{$chr}[$j]->{END} < $leftbound){next};
            if ($self->{NONCONSERVED}{$chr}[$j]->{START} > $rightbound){next};  
            if ($self->{NONCONSERVED}{$chr}[$j]->{START} < $leftbound){
                $maskstart = $leftbound;
            };
            if ($self->{NONCONSERVED}{$chr}[$j]->{END} > $rightbound){
                $maskend = $rightbound; 
            };
            my $element = {START => $maskstart, END => $maskend};
            push @{$self->{GENOME}{$gene[$i]}->{NONCONSERVED}}, $element;
        };
        if ($count > 1000){print'.';$count = 0};
        $count++;
    };
    print"\n\n";
};

sub find_conserved{
    my $self = shift;
    my @gene = keys %{$self->{GENOME}};
    print"Finding conserved regions\n";
    my $count = 0;
    for my $i (0 .. $#gene){
        my ($leftbound, $rightbound) = undef;
        if ($self->{GENOME}{$gene[$i]}{STRAND} =~ /\+/){
            $leftbound = $self->{GENOME}{$gene[$i]}->{TXSTART} - $self->{INFO}{UPSTREAM};
            $rightbound = $self->{GENOME}{$gene[$i]}->{TXEND} + $self->{INFO}{UPSTREAM} - 1;
        }
        else{
            $leftbound = $self->{GENOME}{$gene[$i]}->{TXSTART} - $self->{INFO}{UPSTREAM} + 1;
            $rightbound = $self->{GENOME}{$gene[$i]}->{TXEND} + $self->{INFO}{UPSTREAM};
        };
        my $chr = $self->{GENOME}{$gene[$i]}->{CHROMOSOME};
        my ($up, $down) = $self->find_conserved_bounds($chr, $leftbound, $rightbound);
        for my $j ($up .. $down){
            my $maskstart = $self->{CONSERVED}{$chr}[$j]->{START};
            my $maskend = $self->{CONSERVED}{$chr}[$j]->{END};
            if ($self->{CONSERVED}{$chr}[$j]->{END} < $leftbound){next};
            if ($self->{CONSERVED}{$chr}[$j]->{START} > $rightbound){next};
            if ($self->{CONSERVED}{$chr}[$j]->{START} < $leftbound){
                $maskstart = $leftbound;
            };
            if ($self->{CONSERVED}{$chr}[$j]->{END} > $rightbound){
                $maskend = $rightbound; 
            };
            my $element = {START => $maskstart, END => $maskend};
            push @{$self->{GENOME}{$gene[$i]}->{CONSERVED}}, $element;
        };
        if ($count > 1000){print'.';$count = 0};
        $count++;
    };
    print"\n\n";
};

sub end_chromosome{
    my $self = shift;
    my $chr = shift;
    my $end = 0;
    my $chromfile = "$self->{DIRECTORIES}{GENOME}\/chr$chr\.fa";
    my @args = stat($chromfile);
    my $bytes = $args[7];
    $/ = "\n";
    open(SEQ, "< $chromfile");
    $. = 0; #finds the length of the first line
    my $first_line = <SEQ>;
    my $length_first_line = length($first_line);
    $. = 1; #finds the length of the second line and the record separator
    my $second_line = <SEQ>;
    my $length_full_second_line = length($second_line);
    chomp $second_line;
    my $length_second_line = length ($second_line);
    close SEQ;
    $end = int(($bytes - $length_first_line) * $length_second_line /$length_full_second_line);
    return $end;
};

sub load_refgene{
    my $self = shift;
    my $name = shift;
    my $file = "$self->{DIRECTORIES}{GENOME}\/refGene\.txt";
    print "Finding overlaps in\:\n$file\n";
    my %Records = ();
    open (FH, "< $file");
    while (<FH>){
        my @rec = split(/\t/, $_);
        $rec[1] =~ s/chr//;
        if ($rec[1] !~ /[\dXY]+/){next};
        if ($rec[1] =~ /random/){next};
        if (exists $self->{GENOME}{$rec[0]}{COUNT}){
            my @exonstart = split(/,/, $rec[8]);
            chomp $rec[9];
            my @exonend = split(/,/, $rec[9]);
            my $position = {NAME => $rec[0],
                            CHROMOSOME => $rec[1],
                            STRAND => $rec[2],
                            TXSTART => $rec[3],
                            TXEND => $rec[4],
                            CDSSTART => $rec[5],
                            CDSEND => $rec[6],
                            EXONSTART => \@exonstart,
                            EXONEND => \@exonend,
                            OVERLAPS => undef
                            };
            push @{$Records{$rec[1]}}, $position;
        };
    };
    close FH;
    
    print "Sorting records by transcriptional start sites\n";
    my @Chromosomes = keys %Records;
    @Chromosomes = sort @Chromosomes;
    $self->{CHROMOSOMES} = \@Chromosomes;
    foreach my $Chromosome(@Chromosomes){
        @{$Records{$Chromosome}} = sort { $a->{TXSTART} <=> $b->{TXSTART} } @{$Records{$Chromosome}};
    };
    
    print "Assigning upstream distances\n";
    my $count = 0;
    foreach my $Chromosome(@Chromosomes){
        GENES: for ( my $i = 0; $i <= $#{$Records{$Chromosome}}; $i++ ){
            my ($leftbound, $rightbound) = undef;
            if ($Records{$Chromosome}[$i]->{STRAND} =~ /\+/){
                $leftbound = $Records{$Chromosome}[$i]->{TXSTART} - $self->{INFO}{UPSTREAM};
                $rightbound = $Records{$Chromosome}[$i]->{TXSTART} + $self->{INFO}{DOWNSTREAM} - 1;
            }
            else{
                $leftbound = $Records{$Chromosome}[$i]->{TXEND} - $self->{INFO}{DOWNSTREAM} + 1;
                $rightbound = $Records{$Chromosome}[$i]->{TXEND} + $self->{INFO}{UPSTREAM};
            };
            my $up = $i - 10;
            my $down = $i + 10;
            if ($up < 0){$up = 0};
            if ($down > $#{$Records{$Chromosome}}){$down = $#{$Records{$Chromosome}}};
            EXONS: for (my $j = $up; $j < $down; $j++){
                if ($j == $i){next EXONS};
                if ($Records{$Chromosome}[$j]->{TXEND} < $leftbound){next EXONS};
                if ($Records{$Chromosome}[$j]->{TXSTART} > $rightbound){next EXONS};
                for (my $k = 0; $k <= $#{$Records{$Chromosome}[$j]->{EXONSTART}}; $k++){
                    if ($Records{$Chromosome}[$j]->{EXONEND}[$k] < $leftbound){next EXONS};
                    if ($Records{$Chromosome}[$j]->{EXONSTART}[$k] > $rightbound){next EXONS};
                    my $start = $Records{$Chromosome}[$j]->{EXONSTART}[$k];
                    my $end = $Records{$Chromosome}[$j]->{EXONEND}[$k];
                    if ($Records{$Chromosome}[$j]->{EXONSTART}[$k] < $leftbound){$start = $leftbound};
                    if ($Records{$Chromosome}[$j]->{EXONEND}[$k] > $rightbound){$end = $rightbound};
                    my $element = {START => $start, END => $end, STRAND => $Records{$Chromosome}[$j]->{STRAND}};
                    my $gene = $Records{$Chromosome}[$j]->{NAME};
                    push @{$Records{$Chromosome}[$i]->{OVERLAPS}{$gene}}, $element;
                };
            };
            $count++;
            if ($count > 1000){print'.';$count = 0};
        };
    };
    print"\n";
 
    print "Converting gene records to hash keyed on gene name\n\n";
    my %hash = ();
    my @chroms = keys %Records;
    for my $i (0 .. $#chroms){
        for my $j (0 .. $#{$Records{$chroms[$i]}}){
            $hash{$Records{$chroms[$i]}[$j]->{NAME}} = {
                CHROMOSOME => $Records{$chroms[$i]}[$j]->{CHROMOSOME},
                TXSTART => $Records{$chroms[$i]}[$j]->{TXSTART},
                TXEND => $Records{$chroms[$i]}[$j]->{TXEND},
                STRAND => $Records{$chroms[$i]}[$j]->{STRAND},
                CDSSTART => $Records{$chroms[$i]}[$j]->{CDSSTART},
                CDSEND => $Records{$chroms[$i]}[$j]->{CDSEND},
                EXONSTART => $Records{$chroms[$i]}[$j]->{EXONSTART},
                EXONEND => $Records{$chroms[$i]}[$j]->{EXONEND},
                OVERLAPS => $Records{$chroms[$i]}[$j]->{OVERLAPS},
                CONSERVED => undef,
                NONCONSERVED => undef,
                NONREPEAT => undef,
                KNOWNTFS => undef,
                PRIMERS => undef,
                SCANACE => undef
            };
        };  
    };
    $self->{GENOME} = \%hash;
};

sub load_factors{
    my $self = shift;
    my $organism = shift;
    my $file = shift;
    print "Finding transcription factor sites in\:\n$file\n";
    my %Records = ();
    open (FH, "< $file");
    while (<FH>){
        my @rec = split(/\t/, $_);
        if ($rec[7] !~ /$organism/){next};
        $rec[1] =~ s/chr//;
        if ($rec[1] !~ /[\dXY]+/){next};
        if ($rec[1] =~ /random/){next};
        my $position = {START => $rec[2],
                        END => $rec[3],
                        };
        push @{$Records{$rec[1]}{$rec[8]}}, $position;
    };
    close FH;
    $self->{FACTORS} = \%Records;
    return;
};

sub load_factors_short{
    my $self = shift;
    my $file = shift;
    print "Finding transcription factor sites in\:\n$file\n";
    my %Records = ();
    my %list = ();
    open (FH, "< $file");
    while (<FH>){
        my @rec = split(/\t/, $_);
        $rec[3] =~ s/chr//;
        if ($rec[3] !~ /[\dXY]+/){next};
        if ($rec[3] =~ /random/){next};
        chomp $rec[5];
        my $position = {START => $rec[4],
                        END => $rec[5],
                        FACTOR => $rec[1]
                        };
        push @{$Records{$rec[3]}}, $position;
        $list{$rec[1]}++;
    };
    close FH;
    foreach my $chrom (keys %Records){
        @{$Records{$chrom}} = sort { $a->{START} <=> $b->{START} } @{$Records{$chrom}};
    };
    $Records{LIST} = \%list;
    $self->{FACTORS} = \%Records;
    return;
};

sub load_refgene_all{
    my $self = shift;
    my $name = shift;
    my $file = "$self->{DIRECTORIES}{GENOME}\/refGene\.txt";
    print "Finding overlaps in\:\n$file\n";
    my %Records = ();
    open (FH, "< $file");
    while (<FH>){
        my @rec = split(/\t/, $_);
        $rec[1] =~ s/chr//;
        if ($rec[1] !~ /[\dXY]+/){next};
        if ($rec[1] =~ /random/){next};
        my @exonstart = split(/,/, $rec[8]);
        chomp $rec[9];
        my @exonend = split(/,/, $rec[9]);
        my $position = {NAME => $rec[0],
                        CHROMOSOME => $rec[1],
                        STRAND => $rec[2],
                        TXSTART => $rec[3],
                        TXEND => $rec[4],
                        CDSSTART => $rec[5],
                        CDSEND => $rec[6],
                        EXONSTART => \@exonstart,
                        EXONEND => \@exonend,
                        OVERLAPS => undef
                        };
        push @{$Records{$rec[1]}}, $position;
    };
    close FH;
    
    print "Sorting records by transcriptional start sites\n";
    my @Chromosomes = keys %Records;
    @Chromosomes = sort @Chromosomes;
    $self->{CHROMOSOMES} = \@Chromosomes;
    foreach my $Chromosome(@Chromosomes){
        @{$Records{$Chromosome}} = sort { $a->{TXSTART} <=> $b->{TXSTART} } @{$Records{$Chromosome}};
    };
    
    print "Assigning upstream distances\n";
    my $count = 0;
    foreach my $Chromosome(@Chromosomes){
        GENES: for ( my $i = 0; $i <= $#{$Records{$Chromosome}}; $i++ ){
            my ($leftbound, $rightbound) = undef;
            if ($Records{$Chromosome}[$i]->{STRAND} =~ /\+/){
                $leftbound = $Records{$Chromosome}[$i]->{TXSTART} - $self->{INFO}{UPSTREAM};
                $rightbound = $Records{$Chromosome}[$i]->{TXSTART} + $self->{INFO}{DOWNSTREAM} - 1;
            }
            else{
                $leftbound = $Records{$Chromosome}[$i]->{TXEND} - $self->{INFO}{DOWNSTREAM} + 1;
                $rightbound = $Records{$Chromosome}[$i]->{TXEND} + $self->{INFO}{UPSTREAM};
            };
            my $up = $i - 10;
            my $down = $i + 10;
            if ($up < 0){$up = 0};
            if ($down > $#{$Records{$Chromosome}}){$down = $#{$Records{$Chromosome}}};
            EXONS: for (my $j = $up; $j < $down; $j++){
                if ($j == $i){next EXONS};
                if ($Records{$Chromosome}[$j]->{TXEND} < $leftbound){next EXONS};
                if ($Records{$Chromosome}[$j]->{TXSTART} > $rightbound){next EXONS};
                for (my $k = 0; $k <= $#{$Records{$Chromosome}[$j]->{EXONSTART}}; $k++){
                    if ($Records{$Chromosome}[$j]->{EXONEND}[$k] < $leftbound){next EXONS};
                    if ($Records{$Chromosome}[$j]->{EXONSTART}[$k] > $rightbound){next EXONS};
                    my $start = $Records{$Chromosome}[$j]->{EXONSTART}[$k];
                    my $end = $Records{$Chromosome}[$j]->{EXONEND}[$k];
                    if ($Records{$Chromosome}[$j]->{EXONSTART}[$k] < $leftbound){$start = $leftbound};
                    if ($Records{$Chromosome}[$j]->{EXONEND}[$k] > $rightbound){$end = $rightbound};
                    my $element = {START => $start, END => $end, STRAND => $Records{$Chromosome}[$j]->{STRAND}};
                    my $gene = $Records{$Chromosome}[$j]->{NAME};
                    push @{$Records{$Chromosome}[$i]->{OVERLAPS}{$gene}}, $element;
                };
            };
            $count++;
            if ($count > 1000){print'.';$count = 0};
        };
    };
    print"\n";
 
    print "Converting gene records to hash keyed on gene name\n\n";
    my %hash = ();
    my @chroms = keys %Records;
    for my $i (0 .. $#chroms){
        for my $j (0 .. $#{$Records{$chroms[$i]}}){
            $hash{$Records{$chroms[$i]}[$j]->{NAME}} = {
                CHROMOSOME => $Records{$chroms[$i]}[$j]->{CHROMOSOME},
                TXSTART => $Records{$chroms[$i]}[$j]->{TXSTART},
                TXEND => $Records{$chroms[$i]}[$j]->{TXEND},
                STRAND => $Records{$chroms[$i]}[$j]->{STRAND},
                CDSSTART => $Records{$chroms[$i]}[$j]->{CDSSTART},
                CDSEND => $Records{$chroms[$i]}[$j]->{CDSEND},
                EXONSTART => $Records{$chroms[$i]}[$j]->{EXONSTART},
                EXONEND => $Records{$chroms[$i]}[$j]->{EXONEND},
                OVERLAPS => $Records{$chroms[$i]}[$j]->{OVERLAPS},
                CONSERVED => undef,
                NONCONSERVED => undef,
                NONREPEAT => undef,
                KNOWNTFS => undef,
                PRIMERS => undef,
                SCANACE => undef
            };
        };  
    };
    $self->{GENOME} = \%hash;
};

sub load_conserved{
    my $self = shift;
    my $count = 0;
    for my $i (0 .. $#{$self->{CHROMOSOMES}}){
        my $file = "$self->{DIRECTORIES}{GENOME}\/chr$self->{CHROMOSOMES}[$i]\.axt";
        open (INFILE, "< $file");
        print "Loading conserved regions from\n$file\n";
        my @array = ();
        while (<INFILE>)
        {
            if ($_ =~ /^\d+\s+chr\S+\s+(\d+)\s+(\d+)/)
            {
                my $start = $1;
                my $end = $2;
                my $element = {START => $start, END => $end};
                push @array, $element;
                if ($count > 10000){print'.';$count = 0};
                $count++;
            };
        };
        close INFILE;
        $self->{CONSERVED}{$self->{CHROMOSOMES}[$i]} = \@array;
        print "\n\n";
    };
    my @Chromosomes = @{$self->{CHROMOSOMES}};
    foreach my $Chromosome(@Chromosomes){
        @{$self->{CONSERVED}{$Chromosome}} = sort { $a->{START} <=> $b->{START} } @{$self->{CONSERVED}{$Chromosome}};
        print "Sorting chromosome\: $Chromosome\n";
        my @holder = ();
        my $first = {START => '0',
                    END => $self->{CONSERVED}{$Chromosome}[0]->{START} - 1};
        push @holder, $first;
        for my $j (0 .. $#{$self->{CONSERVED}{$Chromosome}} - 1){
            my $element = {START => ($self->{CONSERVED}{$Chromosome}[$j]->{END} + 1 ),
                           END => ($self->{CONSERVED}{$Chromosome}[$j + 1]->{START} -1 )};
            push @holder, $element;
        }
        my $last = {START => ($self->{CONSERVED}{$Chromosome}[$#{$self->{CONSERVED}{$Chromosome}}]->{END} + 1),
                    END => $self->end_chromosome($Chromosome)};
        push @holder, $last;
        $self->{NONCONSERVED}{$Chromosome} = \@holder;
    };
};

sub load_gene_set{
    my ($self, $name, $infile) = @_;
    my %hash = ();
    print "Loading records\:\n$infile\n";
    open (INFILE, "< $infile");
    while (<INFILE>){
        if ($_ =~ /^(\S+)/){
            my $gene = $1;
            if (exists $self->{GENOME}{$gene}){
                print "$gene\tfound\n";
                $hash{$gene} = $self->{GENOME}{$gene}
            }
            else{print "$gene\tnot found\n"};
        }
    };
    close INFILE;
    print "\n\n";
    $self->{SETS}{$name} = \%hash;
};

sub load_search_set{
    my ($self, $infile) = @_;
    my %hash = ();
    print "Loading records\:\n$infile\n";
    open (INFILE, "< $infile");
    while (<INFILE>){
        if ($_ =~ /^(\S+)/){
            my $gene = $1;
            if (exists $self->{GENOME}{$gene}){
                print "$gene\tfound\n";
                $hash{$gene} = $self->{GENOME}{$gene}
            }
            else{print "$gene\tnot found\n"};
        }
    };
    close INFILE;
    print "\n\n";
    $self->{SEARCHSET} = \%hash;
};

sub load_unique_unigene {
    my $self = shift;
    my $infile = "$self->{DIRECTORIES}{GENOME}\/UniqueUnigene\.txt";
    my %hash = ();
    my $count = 0;
    print "Loading unique UniGene records\:\n$infile\n";
    open (INFILE, "< $infile");
    while (<INFILE>)
    {
        if ($_ =~ /^\>/){
            if ($_ =~ /(NM\_\d+)/){
                $hash{$1}{COUNT}++;
                $count++;
            };
        };
        if ($count > 1000){print'.';$count = 0};
    };
    close INFILE;
    print "\n\n";
    $self->{GENOME} = \%hash;
    my @genes = keys %{$self->{GENOME}};
    $self->{GENES} = \@genes;
};

sub load_unique_unigene_with_testset {
    my $self = shift;
    my $testset = shift;
    my $infile = "$self->{DIRECTORIES}{GENOME}\/UniqueUnigene\.txt";
    my %hash = ();
    my $count = 0;
    print "Loading unique UniGene records\:\n$infile\n";
    open (INFILE, "< $infile");
    while (<INFILE>)
    {
        if ($_ =~ /^\>/){
            if ($_ =~ /(NM\_\d+)/){
                $hash{$1}{COUNT}++;
                $count++;
            };
        };
        if ($count > 1000){print'.';$count = 0};
    };
    close INFILE;
    print "\n\n";
    
    print "Loading records\:\n$testset\n";
    open (INFILE, "< $testset");
    while (<INFILE>){
        if ($_ =~ /(NM\_\d+)/){
            $hash{$1}{COUNT}++;
        }
    };
    close INFILE;
    print "\n\n";
    
    $self->{GENOME} = \%hash;
    my @genes = keys %{$self->{GENOME}};
    $self->{GENES} = \@genes;
};

sub get_scanace_sites {
    my ($self, $set ) = @_;
    my $motifs = $self->{MOTIFS}{$set}->get_number_motifs();
    for my $motif (0 .. $motifs - 1){
        my %sites = $self->{MOTIFS}{$set}->get_motif_sites($motif);
        my @gene = keys %{$self->{SETS}{$set}};
        for my $i (0 .. $#gene){
            for my $j (0 .. $#{$sites{$gene[$i]}{ELEMENT}}){
                my $relstart = $sites{$gene[$i]}{ELEMENT}[$j]->{POSITION};
                my $relend = undef;
                if (!defined $self->{SETS}{$set}{$gene[$i]}{STRAND}){next};
                if($self->{SETS}{$set}{$gene[$i]}{STRAND} =~ /\+/){
                $relstart = $self->{SETS}{$set}{$gene[$i]}{TXSTART} - $self->{INFO}{UPSTREAM} + $relstart;
                    $relend = $relstart + length($sites{$gene[$i]}{ELEMENT}[$j]->{SEQUENCE});
                }
                else{
                    $relstart = $self->{SETS}{$set}{$gene[$i]}{TXEND} + $self->{INFO}{UPSTREAM} - $relstart - length($sites{$gene[$i]}{ELEMENT}[$j]->{SEQUENCE});
                    $relend = $relstart + length($sites{$gene[$i]}{ELEMENT}[$j]->{SEQUENCE});
                };
                my $element = {START => $relstart,
                                STRAND => $sites{$gene[$i]}{ELEMENT}[$j]->{STRAND},
                                SCORE => $sites{$gene[$i]}{ELEMENT}[$j]->{SCORE},
                                END => $relend};
                push @{$self->{SETS}{$set}{$gene[$i]}->{SCANACE}{$motif}}, $element;
            };
        };
    };
};

sub get_factor_sites {
    my ($self, $set ) = @_;
    my @gene = keys %{$self->{SETS}{$set}};
    for my $i (0 .. $#gene){
        my ($leftbound, $rightbound) = undef;
        if ($self->{SETS}{$set}{$gene[$i]}{STRAND} =~ /\+/){
            $leftbound = $self->{SETS}{$set}{$gene[$i]}->{TXSTART} - $self->{INFO}{UPSTREAM};
            $rightbound = $self->{SETS}{$set}{$gene[$i]}->{TXSTART} + $self->{INFO}{DOWNSTREAM} - 1;
        }
        else{
            $leftbound = $self->{SETS}{$set}{$gene[$i]}->{TXEND} - $self->{INFO}{DOWNSTREAM} + 1;
            $rightbound = $self->{SETS}{$set}{$gene[$i]}->{TXEND} + $self->{INFO}{UPSTREAM};
        };
        my $chr = $self->{SETS}{$set}{$gene[$i]}->{CHROMOSOME};
        my ($startarray, $endarray) = $self->find_factor_bounds($chr, $leftbound, $rightbound);
        for my $j ($startarray .. $endarray){
            my $element = {START => $self->{FACTORS}{$chr}[$j]->{START},
                            END => $self->{FACTORS}{$chr}[$j]->{END}};
            push @{$self->{SETS}{$set}{$gene[$i]}->{FACTOR}{$self->{FACTORS}{$chr}[$j]->{FACTOR}}}, $element;
        };
    };
};

sub make_randoms{
    my ($self, $file, $sets) = @_;
    my %genes = ();
    $self->read_sequences($file);
    $self->sort_sequences();
    $file =~ s/\.\w+$//;
    for my $i (1 .. $sets){
        my $out = "$file"."random$i\.seq";
        print"Making random set\:\n$out\n";
        open (OUT, "> $out");
        for my $j (0 .. $#{$self->{SEQUENCES}}){
            my ($gene, $seq) = $self->get_randoms($self->{SEQUENCES}[$j]->{LENGTH}, \%genes);
            $genes{$gene}++;
            print OUT "\>$gene\n$seq\n";
        };
        close OUT;
    };
};

sub make_upstream_randoms {
    my ($self, $file, $sets) = @_;
    my %genes = ();
    $self->read_sequences($file);
    $self->sort_sequences();
    $file =~ s/\.\w+$//;
    for my $i (1 .. $sets){
        my $out = "$file"."$i"."Random\.seq";
        print"Making random set\:\n$out\n";
        open (OUT, "> $out");
        for my $j (0 .. $#{$self->{SEQUENCES}}){
            my ($gene, $seq) = $self->get_upstream_randoms($self->{SEQUENCES}[$j]->{LENGTH}, \%genes);
            $genes{$gene}++;
            print OUT "\>$gene\_1\n$seq\n";
        };
        close OUT;
    };
};

sub make_first_intron_randoms {
    my ($self, $file, $sets) = @_;
    my %genes = ();
    $self->read_sequences($file);
    $self->sort_sequences();
    $file =~ s/\.\w+$//;
    for my $i (1 .. $sets){
        my $out = "$file"."$i"."Random\.seq";
        print"Making random set\:\n$out\n";
        open (OUT, "> $out");
        for my $j (0 .. $#{$self->{SEQUENCES}}){
            my ($gene, $seq) = $self->get_first_intron_randoms($self->{SEQUENCES}[$j]->{LENGTH}, \%genes);
            $genes{$gene}++;
            print OUT "\>$gene\_1\n$seq\n";
        };
        close OUT;
    };
};

sub make_downstream_randoms {
    my ($self, $file, $sets) = @_;
    my %genes = ();
    $self->read_sequences($file);
    $self->sort_sequences();
    $file =~ s/\.\w+$//;
    for my $i (1 .. $sets){
        my $out = "$file"."$i"."Random\.seq";
        print"Making random set\:\n$out\n";
        open (OUT, "> $out");
        for my $j (0 .. $#{$self->{SEQUENCES}}){
            my ($gene, $seq) = $self->get_downstream_randoms($self->{SEQUENCES}[$j]->{LENGTH}, \%genes);
            $genes{$gene}++;
            print OUT "\>$gene\_1\n$seq\n";
        };
        close OUT;
    };
};

sub sort_sequences{
    my $self = shift;
    @{$self->{SEQUENCES}} = sort { $b->{LENGTH} <=> $a->{LENGTH}} @{$self->{SEQUENCES}};
};

sub read_sequences{
    my ($self, $file) = @_;
    $self->{SEQUENCES} = ();
    open (FILE, "< $file");
    $/ = '>';
    my @array = ();
    while (<FILE>){push @array, [split(/\n/, $_)]};
    close FILE;
    $/ = "\n";
    for my $i (1 .. $#array){
        if (length($array[$i][0]) =~ /\>/){next};
        if (length($array[$i][1]) < 1){next};
        my $element = { NAME => $array[$i][0], SEQUENCE => $array[$i][1], LENGTH => length($array[$i][1])};
        push @{$self->{SEQUENCES}}, $element;
    };
};

sub get_randoms {
    my ($self, $size, $hashref) = @_;    
    my $test = '1';
    my $seq = '';
    my $gene = '';
    my %genes = %{$hashref};
    while ($test){
        $gene = int(rand($#{$self->{GENES}}));
        if (exists $genes{$self->{GENES}[$gene]}){next};
        my ($leftbound, $rightbound) = undef;
        if ($self->{GENOME}{$self->{GENES}[$gene]}{STRAND} =~ /\+/){
            $leftbound = $self->{GENOME}{$self->{GENES}[$gene]}->{TXSTART} - $self->{INFO}{UPSTREAM};
            $rightbound = $self->{GENOME}{$self->{GENES}[$gene]}->{TXSTART} + $self->{INFO}{DOWNSTREAM};
        }
        else{
            $leftbound = $self->{GENOME}{$self->{GENES}[$gene]}->{TXEND} - $self->{INFO}{DOWNSTREAM};
            $rightbound = $self->{GENOME}{$self->{GENES}[$gene]}->{TXEND} + $self->{INFO}{UPSTREAM};
        };
        my $chr = $self->{GENOME}{$self->{GENES}[$gene]}->{CHROMOSOME};
        if ($chr !~ /\S+/){next};
        my $sequence = $self->get_sequence($chr, $leftbound, $rightbound);
        $sequence = $self->mask_overlaps($chr, $leftbound, $rightbound, $self->{GENES}[$gene], $sequence);
        if ($self->{INFO}{CONSERVED}){
            $sequence = $self->mask_nonconserved($self->{GENOME}{$self->{GENES}[$gene]}->{CHROMOSOME}, $leftbound, $rightbound, $sequence);
        };
        if ($self->{INFO}{MASKED}){$sequence = mask_repeats($sequence)};
        $self->{GENOME}{$self->{GENES}[$gene]}->{SEQUENCES} = split_sequence($sequence, $leftbound);
        for my $j (0 .. $#{$self->{GENOME}{$self->{GENES}[$gene]}->{SEQUENCES}}){
            $seq = $self->{GENOME}{$self->{GENES}[$gene]}->{SEQUENCES}[$j]->{SEQUENCE};
            $seq = uppercase_sequence($seq);
            if (length($seq) < $size){next};
            if (length($seq) > $size){$seq = substr($seq, 0, $size)};
            if (length($seq) == $size){$test = 0; last};
        };
    };
    return ($self->{GENES}[$gene], $seq);
};

sub get_upstream_randoms {
    my ($self, $size, $hashref) = @_;    
    my $test = '1';
    my $seq = '';
    my $gene = '';
    my %genes = %{$hashref};
    while ($test){
        $gene = int(rand($#{$self->{GENES}}));
        if (exists $genes{$self->{GENES}[$gene]}){next};
        my ($leftbound, $rightbound) = undef;
        if (!defined $self->{GENOME}{$self->{GENES}[$gene]}{STRAND}){next};
        if ($self->{GENOME}{$self->{GENES}[$gene]}{STRAND} =~ /\+/){
            $leftbound = $self->{GENOME}{$self->{GENES}[$gene]}->{TXSTART} - $self->{INFO}{UPSTREAM};
            $rightbound = $self->{GENOME}{$self->{GENES}[$gene]}->{TXSTART} + $self->{INFO}{DOWNSTREAM} - 1;
        }
        else{
            $leftbound = $self->{GENOME}{$self->{GENES}[$gene]}->{TXEND} - $self->{INFO}{DOWNSTREAM} + 1;
            $rightbound = $self->{GENOME}{$self->{GENES}[$gene]}->{TXEND} + $self->{INFO}{UPSTREAM};
        };
        my $chr = $self->{GENOME}{$self->{GENES}[$gene]}->{CHROMOSOME};
        if ($chr !~ /\S+/){next};
        my $sequence = $self->get_sequence($chr, $leftbound, $rightbound);
        $sequence = $self->mask_overlaps($chr, $leftbound, $rightbound, $self->{GENES}[$gene], $sequence);
        if ($self->{INFO}{CONSERVED}){
            $sequence = $self->mask_nonconserved($self->{GENOME}{$self->{GENES}[$gene]}->{CHROMOSOME}, $leftbound, $rightbound, $sequence);
        };
        if ($self->{INFO}{MASKED}){$sequence = mask_repeats($sequence)};
        my @seqs = @{split_sequence($sequence, $leftbound)};
        my $j =  int(rand($#seqs + 1));
        if (!defined $seqs[$j]->{SEQUENCE}){next}; 
        if ($self->{GENOME}{$self->{GENES}[$gene]}->{STRAND} =~ /\+/){
            $seq = $seqs[$j]->{SEQUENCE};
            $seq = uppercase_sequence($seq);
            if (length($seq) < $size){next};
            if (length($seq) > $size){
                my $begin = int(rand(length($seq) - $size));
                $seq = substr($seq, $begin, $size)
                };
            if (length($seq) == $size){$test = 0; last};
        }
        else{
            $seq = $seqs[$j]->{SEQUENCE};
            $seq = reverse_complement($seq);
            $seq = uppercase_sequence($seq);
            if (length($seq) < $size){next};
            if (length($seq) > $size){
                my $begin = int(rand(length($seq) - $size));
                $seq = substr($seq, $begin, $size)
                };
            if (length($seq) == $size){$test = 0; last};
        };
    };
    return ($self->{GENES}[$gene], $seq);
};

sub get_first_intron_randoms {
    my ($self, $size, $hashref) = @_;    
    my $test = '1';
    my $seq = '';
    my $gene = '';
    my %genes = %{$hashref};
    while ($test){
        $gene = int(rand($#{$self->{GENES}}));
        if (exists $genes{$self->{GENES}[$gene]}){next};
        my ($leftbound, $rightbound) = undef;
        
        if (!exists $self->{GENOME}{$self->{GENES}[$gene]}->{EXONSTART}[1]){next};
        if ($self->{GENOME}{$self->{GENES}[$gene]}{STRAND} =~ /\+/){
            $leftbound = $self->{GENOME}{$self->{GENES}[$gene]}->{EXONEND}[0] + 1;
            $rightbound = $self->{GENOME}{$self->{GENES}[$gene]}->{EXONSTART}[1] - 1;
            if (($rightbound - $leftbound) > $self->{INFO}{UPSTREAM}){$rightbound = $leftbound + $self->{INFO}{UPSTREAM}};
        }
        else{
            $leftbound = $self->{GENOME}{$self->{GENES}[$gene]}->{EXONEND}[$#{$self->{GENOME}{$self->{GENES}[$gene]}->{EXONEND}} - 1] + 1;
            $rightbound = $self->{GENOME}{$self->{GENES}[$gene]}->{EXONSTART}[$#{$self->{GENOME}{$self->{GENES}[$gene]}->{EXONSTART}}] - 1;
            if (($rightbound - $leftbound) > $self->{INFO}{UPSTREAM}){$leftbound = $rightbound - $self->{INFO}{UPSTREAM}};
        };
        my $chr = $self->{GENOME}{$self->{GENES}[$gene]}->{CHROMOSOME};
        if ($chr !~ /\S+/){next};
        my $sequence = $self->get_sequence($chr, $leftbound, $rightbound);
        $sequence = $self->mask_overlaps($chr, $leftbound, $rightbound, $self->{GENES}[$gene], $sequence);
        if ($self->{INFO}{CONSERVED}){
            $sequence = $self->mask_nonconserved($self->{GENOME}{$self->{GENES}[$gene]}->{CHROMOSOME}, $leftbound, $rightbound, $sequence);
        };
        if ($self->{INFO}{MASKED}){$sequence = mask_repeats($sequence)};
        my @seqs = @{split_sequence($sequence, $leftbound)};
        my $j =  int(rand($#seqs + 1));
        if (!defined $seqs[$j]->{SEQUENCE}){next};
        if ($self->{GENOME}{$self->{GENES}[$gene]}->{STRAND} =~ /\+/){
            $seq = $seqs[$j]->{SEQUENCE};
            $seq = uppercase_sequence($seq);
            if (length($seq) < $size){next};
            if (length($seq) > $size){
                my $begin = int(rand(length($seq) - $size));
                $seq = substr($seq, $begin, $size)
                };
            if (length($seq) == $size){$test = 0; last};
        }
        else{
            $seq = $seqs[$j]->{SEQUENCE};
            $seq = reverse_complement($seq);
            $seq = uppercase_sequence($seq);
            if (length($seq) < $size){next};
            if (length($seq) > $size){
                my $begin = int(rand(length($seq) - $size));
                $seq = substr($seq, $begin, $size)
                };
            if (length($seq) == $size){$test = 0; last};
        };
    };
    return ($self->{GENES}[$gene], $seq);
};

sub get_downstream_randoms {
    my ($self, $size, $hashref) = @_;    
    my $test = '1';
    my $seq = '';
    my $gene = '';
    my %genes = %{$hashref};
    while ($test){
        $gene = int(rand($#{$self->{GENES}}));
        if (exists $genes{$self->{GENES}[$gene]}){next};
        if (!exists $self->{GENOME}{$self->{GENES}[$gene]}{STRAND}){next};
        my ($leftbound, $rightbound) = undef;
        if ($self->{GENOME}{$self->{GENES}[$gene]}{STRAND} =~ /\+/){
            $leftbound = $self->{GENOME}{$self->{GENES}[$gene]}->{TXEND} - $self->{INFO}{DOWNSTREAM} + 1;
            $rightbound = $self->{GENOME}{$self->{GENES}[$gene]}->{TXEND} + $self->{INFO}{UPSTREAM};
        }
        else{
            $leftbound = $self->{GENOME}{$self->{GENES}[$gene]}->{TXSTART} - $self->{INFO}{UPSTREAM};
            $rightbound = $self->{GENOME}{$self->{GENES}[$gene]}->{TXSTART} + $self->{INFO}{DOWNSTREAM} - 1;
        };
        my $chr = $self->{GENOME}{$self->{GENES}[$gene]}->{CHROMOSOME};
        if ($chr !~ /\S+/){next};
        my $sequence = $self->get_sequence($chr, $leftbound, $rightbound);
        $sequence = $self->mask_overlaps($chr, $leftbound, $rightbound, $self->{GENES}[$gene], $sequence);
        if ($self->{INFO}{CONSERVED}){
            $sequence = $self->mask_nonconserved($self->{GENOME}{$self->{GENES}[$gene]}->{CHROMOSOME}, $leftbound, $rightbound, $sequence);
        };
        if ($self->{INFO}{MASKED}){$sequence = mask_repeats($sequence)};
        my @seqs = @{split_sequence($sequence, $leftbound)};
        my $j =  int(rand($#seqs + 1));
        if (!defined $seqs[$j]->{SEQUENCE}){next};
        if ($self->{GENOME}{$self->{GENES}[$gene]}->{STRAND} =~ /\+/){
            $seq = $seqs[$j]->{SEQUENCE};
            $seq = uppercase_sequence($seq);
            if (length($seq) < $size){next};
            if (length($seq) > $size){
                my $begin = int(rand(length($seq) - $size));
                $seq = substr($seq, $begin, $size)
                };
            if (length($seq) == $size){$test = 0; last};
        }
        else{
            $seq = $seqs[$j]->{SEQUENCE};
            $seq = reverse_complement($seq);
            $seq = uppercase_sequence($seq);
            if (length($seq) < $size){next};
            if (length($seq) > $size){
                my $begin = int(rand(length($seq) - $size));
                $seq = substr($seq, $begin, $size)
                };
            if (length($seq) == $size){$test = 0; last};
        };
    };
    return ($self->{GENES}[$gene], $seq);
};

sub list_intron_counts{
    my $self = shift;
    my $file = shift;
    my $sitelength = shift;
    my @array = ();
    foreach my $site (keys %{$self->{FREQUENCY}{$sitelength}}){
        my $element = {NAME => $site,
                       COUNT => $self->{FREQUENCY}{$sitelength}{$site}->{COUNT},
                       FREQ => $self->{FREQUENCY}{$sitelength}{$site}->{FREQ}
                        };
        push @array, $element;
    };
    @array = sort { $b->{COUNT} <=> $a->{COUNT} } @array;
    open(FILE, "> $file");
    for my $i (0 .. $#array){
        my $rank = $i + 1;
        print FILE "$rank\t$array[$i]->{NAME}\t$array[$i]->{COUNT}\t$array[$i]->{FREQ}\n";
    };
    close FILE;
};

sub calculate_genome_markov {
    my $self = shift;
    my $sitelength = shift;
    my @chroms = keys %{$self->{INTRONIC}};
    foreach my $chroms (keys %{$self->{INTRONIC}}){
        print "Chromosome\: $chroms\n";
        my $count = 0;
        for my $i (0 .. $#{$self->{INTRONIC}{$chroms}}){
            my ($leftbound, $rightbound) = undef;
            $leftbound = $self->{INTRONIC}{$chroms}[$i]->{START};
            $rightbound = $self->{INTRONIC}{$chroms}[$i]->{END};
            my $sequence = $self->get_sequence($chroms, $leftbound, $rightbound);
            if ($self->{INFO}{CONSERVED}){
            $sequence = $self->mask_nonconserved($chroms, $leftbound, $rightbound, $sequence);
            };
            if ($self->{INFO}{MASKED}){$sequence = mask_repeats($sequence)};
            $sequence = uppercase_sequence($sequence);
            my $seqref = split_sequence($sequence, $leftbound);
            for (my $j = 0; $j <= $#{$seqref}; $j++){
                $self->load_markov($sitelength, $seqref->[$j]->{SEQUENCE});
            };
            if ($count >= 1000){print "\."; $count = 0}else{$count++};
        };
        if ($count >= 1000){print "\."; $count = 0}else{$count++};
        print "\n";
    };
    my $total = 0;
    foreach my $site (keys %{$self->{FREQUENCY}{$sitelength}}){
        $total += $self->{FREQUENCY}{$sitelength}{$site}->{COUNT};
    };
    foreach my $site (keys %{$self->{FREQUENCY}{$sitelength}}){
        $self->{FREQUENCY}{$sitelength}{$site}->{FREQ} = $self->{FREQUENCY}{$sitelength}{$site}->{COUNT} / $total;
    };
};

sub calculate_markov {
    my $self = shift;
    my $file = shift;
    my $order = shift;
    $order++;
    print "Calculating Markov model\n";
    print "$file\n";
    print "Counting words\n";
    my $count = 0;
    $self->make_markov($order);
    open (FASTA, "< $file");
    while (<FASTA>){
        if ($_ =~ /^\>/){next};
        chomp;
        my $seq = $_;
        for my $i (1 .. $order){
            $self->load_markov($i, $seq);
        };
        if ($count > 1000){print "."; $count = 0};
        $count++;
    };
    print "\nCalculating word frequencies\n";
    close FASTA;
    for my $j (1 .. $order){
        my $total = 0;
        foreach my $site (keys %{$self->{FREQUENCY}{$j}}){
            $total += $self->{FREQUENCY}{$j}{$site}->{COUNT};
        };
        foreach my $site (keys %{$self->{FREQUENCY}{$j}}){
            $self->{FREQUENCY}{$j}{$site}->{FREQ} = $self->{FREQUENCY}{$j}{$site}->{COUNT} / $total;
        };
    };
    print"Done\n\n"; 
    return;
};

sub list_markov{
    my $self = shift;
    my $file = shift;
    my $order = shift;
    $order++;
    print "Listing Markov model\n";
    print "$file\n";
    
    open (MARKOV, "> $file");
    for my $i (1 .. $order){
        printf MARKOV ("# order %d\n", $i-1);
        foreach my $tuple (keys %{$self->{FREQUENCY}{$i}}){
            printf MARKOV ("%s %9.3e\n", $tuple, $self->{FREQUENCY}{$i}{$tuple}->{FREQ})
        };
        printf ("Order %d\n", $i-1);
    }
    close MARKOV;
    print "Done\n\n";
    return;
};

sub calculate_counts {
    my $self = shift;
    my $file = shift;
    my $sitelength = shift;
    print "Calculating Counts\n";
    print "$file\n";
    print "Counting words\n";
    my $count = 0;
    $self->make_count_holder($sitelength);
    open (FASTA, "< $file");
    while (<FASTA>){
        if ($_ =~ /^\>/){next};
        chomp;
        my $seq = $_;
        $self->load_markov($sitelength, $seq);
        if ($count > 1000){print "."; $count = 0};
        $count++;
    };
    print "\nCalculating word frequencies\n";
    close FASTA;

    my $total = 0;
    foreach my $site (keys %{$self->{FREQUENCY}{$sitelength}}){
        $total += $self->{FREQUENCY}{$sitelength}{$site}->{COUNT};
    };
    foreach my $site (keys %{$self->{FREQUENCY}{$sitelength}}){
        $self->{FREQUENCY}{$sitelength}{$site}->{FREQ} = $self->{FREQUENCY}{$sitelength}{$site}->{COUNT} / $total;
    };
    print"Done\n\n"; 
    return;
};

sub subtract_counts {
    my $self = shift;
    my $words = shift;
    my $sitelength = shift;
    print "Subtracting Counts\n";
    foreach my $site (keys %{$self->{FREQUENCY}{$sitelength}}){
        $self->{FREQUENCY}{$sitelength}{$site}->{COUNT} = $self->{FREQUENCY}{$sitelength}{$site}->{COUNT} - $words->{FREQUENCY}{$sitelength}{$site}->{COUNT};
    };
    my $total = 0;
    foreach my $site (keys %{$self->{FREQUENCY}{$sitelength}}){
        $total += $self->{FREQUENCY}{$sitelength}{$site}->{COUNT};
    };
    foreach my $site (keys %{$self->{FREQUENCY}{$sitelength}}){
        $self->{FREQUENCY}{$sitelength}{$site}->{FREQ} = $self->{FREQUENCY}{$sitelength}{$site}->{COUNT} / $total;
    };
    print"Done\n\n"; 
    return;
};

sub load_counts{
    my $self = shift;
    my $file = shift;
    open(FILE, "< $file");
    while (<FILE>){
        my ($name, $count, $frequency);
        if ($_ =~ /\d+\s+(\S+)\s+(\d+)\s+(\S+)/){$name = $1; $count = $2; $frequency = $3}
        else{next};
        $self->{FREQUENCY}{length($name)}{$name}->{COUNT} = $count;
        $self->{FREQUENCY}{length($name)}{$name}->{FREQ} = $frequency;
    };
    close FILE;
};

sub list_counts{
    my $self = shift;
    my $file = shift;
    my $sitelength = shift;
    my @array = ();
    foreach my $site (keys %{$self->{FREQUENCY}{$sitelength}}){
        my $element = {NAME => $site,
                       COUNT => $self->{FREQUENCY}{$sitelength}{$site}->{COUNT},
                       FREQ => $self->{FREQUENCY}{$sitelength}{$site}->{FREQ}
                        };
        push @array, $element;
    };
    @array = sort { $b->{COUNT} <=> $a->{COUNT} } @array;
    open(FILE, "> $file");
    for my $i (0 .. $#array){
        my $rank = $i + 1;
        print FILE "$rank\t$array[$i]->{NAME}\t$array[$i]->{COUNT}\t$array[$i]->{FREQ}\n";
    };
    close FILE;
};

sub compare_counts{
    my $self = shift;
    my $comp = shift;
    my $file = shift;
    my $sitelength = shift;
    my @array = ();
    foreach my $site (keys %{$self->{FREQUENCY}{$sitelength}}){
        my $ratio = $self->{FREQUENCY}{$sitelength}{$site}->{FREQ} / $comp->{FREQUENCY}{$sitelength}{$site}->{FREQ};
        my $element = {NAME => $site,
                       FREQ1 => $self->{FREQUENCY}{$sitelength}{$site}->{FREQ},
                       FREQ2 => $comp->{FREQUENCY}{$sitelength}{$site}->{FREQ},
                       RATIO => $ratio
                        };
        push @array, $element;
    };
    @array = sort { $b->{RATIO} <=> $a->{RATIO} } @array;
    open(FILE, "> $file");
    for my $i (0 .. $#array){
        my $rank = $i + 1;
        print FILE "$rank\t$array[$i]->{NAME}\t$array[$i]->{FREQ1}\t$array[$i]->{FREQ2}\t$array[$i]->{RATIO}\n";
    }
    close FILE;
};

sub load_markov {
    my $self = shift;
    my $sitelength = shift;
    my $seq = shift;
    my $sites = length($seq) - $sitelength;
    for my $j(0 .. $sites){
        my $holder = substr($seq, $j, $sitelength);
        if (length($holder) == $sitelength and $holder =~ /^[AGTC]{$sitelength,}$/){
            if (defined $self->{FREQUENCY}{$sitelength}{$holder}->{COUNT}){
                $self->{FREQUENCY}{$sitelength}{$holder}->{COUNT} += 1}
            else{print"$holder not found\n"};
            $holder = reverse($holder);
            $holder =~ tr /AGTC/TCAG/;
            if (defined $self->{FREQUENCY}{$sitelength}{$holder}->{COUNT}){
                $self->{FREQUENCY}{$sitelength}{$holder}->{COUNT} += 1
            }
            else{print"$holder not found\n"};
        };
    };
    
};

sub make_markov {
    my $self = shift;
    my $order = shift;
    for my $i (1 .. $order){
        $self->{FREQUENCY}{$i} = get_base_string($i);
    };
};

sub make_count_holder {
    my $self = shift;
    my $sitelength = shift;
    $self->{FREQUENCY}{$sitelength} = get_base_string($sitelength);
};

sub get_base_string{
    my $depth = shift;
    my %hash = ();
    get_bases(1, '', $depth, \%hash);
    return \%hash
};

sub get_bases {
    my $level = shift;
    my $holder = shift;
    my $depth = shift;
    my $hashref = shift;
    my @base = ();
    $base[0] = 'A';
    $base[1] = 'G';
    $base[2] = 'T';
    $base[3] = 'C';
    for my $i (0 .. 3){
        my $dna = "$holder"."$base[$i]";
        if ($level != $depth){
            get_bases(($level + 1), $dna, $depth, $hashref)
        }
        else{my $element = {COUNT => "0", FREQ => "0"}; ${$hashref}{$dna} = $element;};
    };
};

sub upstream_scan_view
{
    my ($self, $set, $region) = @_;    
    my $count = 0;
    my $radius = 1;
    my $image_width = 1200;
    my $number = $self->{MOTIFS}{$set}->get_number_motifs();
    my $image_height = 50 * $number;

    my $scale = 500/($self->{INFO}{UPSTREAM} + $self->{INFO}{DOWNSTREAM});
    my @gene = keys %{$self->{SEARCHSET}};
    @gene = sort { ($a =~ /(\d+)/)[0] <=> ($b =~ /(\d+)/)[0] } @gene;
    
    for my $i (0 .. $#gene){
        my $file = "$self->{DIRECTORIES}{WORK}\/$set$region$gene[$i]\.svg";
        my $io = new IO::File;
        $io->open ("> $file");
        my $writer = new XML::Writer(OUTPUT => $io, DATA_MODE => 1, DATA_INDENT => 1);
        
        $writer->startTag(
            'svg',
            height => $image_height,
            width  => $image_width
        );
        
        $writer->emptyTag(
            'rect',
            height => $image_height,
            width  => $image_width,
            fill   => 'white'
        );
        
        $writer->startTag(
            'g',
            id => 'motif-page',
            transform => 'translate(100,100)'
        );
    
        my $position = 0;
        if ($self->{SETS}{$set}{$gene[$i]}->{STRAND} =~ /\+/){
            for my $m (0 .. $number - 1){
                for my $p (0 .. $#{$self->{SETS}{$set}{$gene[$i]}->{SCANACE}{$m}}){
                    my $rec_start = $self->{SETS}{$set}{$gene[$i]}->{SCANACE}{$m}[$p]->{START} -
                    $self->{SETS}{$set}{$gene[$i]}->{TXSTART} +
                    $self->{INFO}{UPSTREAM};
                
                    my $rec_end = $self->{SETS}{$set}{$gene[$i]}->{SCANACE}{$m}[$p]->{END} -
                    $self->{SETS}{$set}{$gene[$i]}->{TXSTART} +
                    $self->{INFO}{UPSTREAM};
                
                    my $rec_length = $rec_end - $rec_start;
                    my $line_height = $position + 2;
                    my $line_start = $scale * $rec_start;
                    my $line_length = $scale * $rec_length;
                    if ($line_start > 500){next};
                    if ($line_start + $line_length > 500){$line_length = 500 - $line_start};
                
                    $writer->emptyTag(
                        'rect',
                        x => $line_start,
                        y  => $line_height,
                        width => $line_length,
                        height  => '8',
                        style => 'fill: black;stroke-width: 0.1; stroke: black;'
                    );
                };
                my $pos_text = $position + 10;
                $writer->startTag('text',
                    x=> '525',
                    y=> $pos_text,
                    style => 'font-family:Arial-BoldMT;font-size:12;fill:black;text-anchor:front;'
                );
                my $motnum = $m + 1;
                $writer->characters("Motif $motnum");
                $writer->endTag('text');
                
                $writer->startTag(
                'g',
                id => 'motif-position',
                transform => "translate\(600\,$pos_text\)",
                );
                $writer->startTag(
                'g',
                id => 'motif-scale',
                transform => "scale\(0\.3\)",
                );
                $self->{MOTIFS}{$set}->get_motif_svg($io, $writer, ($m));
                $writer->endTag('g');
                $writer->endTag('g');
                
                $position += 30;
            };
            
            my $offset = $position + 6;
            
            $writer->emptyTag(
                'line',
                x1 => '0',
                y1  => $offset,
                x2 => '500',
                y2  => $offset,
                style => 'stroke-width: 4; stroke: black;'
            );
            
            my $pos_text = $position + 10;
            $writer->startTag('text',
                x=> '525',
                y=> $pos_text,
                style => 'font-family:Arial-BoldMT;font-size:12;fill:black;text-anchor:front;'
            );
            $writer->characters("$gene[$i]");
            $writer->endTag('text');
            
            for my $z (0 .. $#{$self->{SETS}{$set}{$gene[$i]}->{NONREPEAT}}){
                my $rec_start = $self->{SETS}{$set}{$gene[$i]}->{NONREPEAT}[$z]->{START} -
                $self->{SETS}{$set}{$gene[$i]}->{TXSTART} +
                $self->{INFO}{UPSTREAM};
                
                my $rec_end = $self->{SETS}{$set}{$gene[$i]}->{NONREPEAT}[$z]->{END} -
                $self->{SETS}{$set}{$gene[$i]}->{TXSTART} +
                $self->{INFO}{UPSTREAM};
                
                my $rec_length = $rec_end - $rec_start;
                my $line_height = $position + 2;
                my $line_start = $scale * $rec_start;
                my $line_length = $scale * $rec_length;
                if ($line_start > 500){next};
                if ($line_start + $line_length > 500){$line_length = 500 - $line_start};
                
                $writer->emptyTag(
                        'rect',
                        x => $line_start,
                        y  => $line_height,
                        width => $line_length,
                        height  => '8',
                        rx => $radius,
                        style => 'fill: black;stroke-width: 1; stroke: black;'
                ); 
            };
            
            for my $z (0 .. $#{$self->{SETS}{$set}{$gene[$i]}->{CONSERVED}}){
                
                my $rec_start = $self->{SETS}{$set}{$gene[$i]}->{CONSERVED}[$z]->{START} -
                $self->{SETS}{$set}{$gene[$i]}->{TXSTART} +
                $self->{INFO}{UPSTREAM};
                
                my $rec_end = $self->{SETS}{$set}{$gene[$i]}->{CONSERVED}[$z]->{END} -
                $self->{SETS}{$set}{$gene[$i]}->{TXSTART} +
                $self->{INFO}{UPSTREAM};
                
                my $rec_length = $rec_end - $rec_start;
                my $line_height = $position + 2;
                my $line_start = $scale * $rec_start;
                my $line_length = $scale * $rec_length;
                if ($line_start > 500){next};
                if ($line_start + $line_length > 500){$line_length = 500 - $line_start};
                    
                $writer->emptyTag(
                    'rect',
                    x => $line_start,
                    y  => $position,
                    width => $line_length,
                    height  => '12',
                    rx => $radius,
                    style => 'fill: red;stroke-width: 1; stroke: black;'
                ); 
            };
            $position = $position + 30;
        }
        else{
            #add scans here
            for my $m (0 .. $number - 1){
                for my $p (0 .. $#{$self->{SETS}{$set}{$gene[$i]}->{SCANACE}{$m}}){
                    my $rec_start = $self->{SETS}{$set}{$gene[$i]}->{TXEND} +
                    $self->{INFO}{UPSTREAM} -
                    $self->{SETS}{$set}{$gene[$i]}->{SCANACE}{$m}[$p]->{END};
                
                    my $rec_end = 
                    $self->{SETS}{$set}{$gene[$i]}->{TXEND} +
                    $self->{INFO}{UPSTREAM} -
                    $self->{SETS}{$set}{$gene[$i]}->{SCANACE}{$m}[$p]->{START};
                
                    my $rec_length = $rec_end - $rec_start;
                    my $line_height = $position + 2;
                    my $line_start = $scale * $rec_start;
                    my $line_length = $scale * $rec_length;
                    if ($line_start > 500){next};
                    if ($line_start + $line_length > 500){$line_length = 500 - $line_start};
                
                    $writer->emptyTag(
                        'rect',
                        x => $line_start,
                        y  => $line_height,
                        width => $line_length,
                        height  => '8',
                        style => 'fill: black;stroke-width: 0.1; stroke: black;'
                    ); 
                };
                my $pos_text = $position + 10;
                $writer->startTag('text',
                    x=> '525',
                    y=> $pos_text,
                    style => 'font-family:Arial-BoldMT;font-size:12;fill:black;text-anchor:front;'
                );
                my $motnum = $m + 1;
                $writer->characters("Motif $motnum");
                $writer->endTag('text');
                
                $writer->startTag(
                'g',
                id => 'motif-position',
                transform => "translate\(600\,$pos_text\)",
                );
                $writer->startTag(
                'g',
                id => 'motif-scale',
                transform => "scale\(0\.3\)",
                );
                $self->{MOTIFS}{$set}->get_motif_svg($io, $writer, ($m));
                $writer->endTag('g');
                $writer->endTag('g');
                
                $position += 30;
            };
            
            my $offset = $position + 6;
            
            $writer->emptyTag(
                'line',
                x1 => '0',
                y1  => $offset,
                x2 => '500',
                y2  => $offset,
                style => 'stroke-width: 4; stroke: black;'
            );
            
            my $pos_text = $position + 10;
            $writer->startTag('text',
                x=> '525',
                y=> $pos_text,
                style => 'font-family:Arial-BoldMT;font-size:12;fill:black;text-anchor:front;'
            );
            $writer->characters("$gene[$i]");
            $writer->endTag('text');
            
            for my $z (0 .. $#{$self->{SETS}{$set}{$gene[$i]}->{NONREPEAT}}){
                my $rec_end =  $self->{SETS}{$set}{$gene[$i]}->{TXEND} +
                $self->{INFO}{UPSTREAM} -
                $self->{SETS}{$set}{$gene[$i]}->{NONREPEAT}[$z]->{START};
                
                my $rec_start = $self->{SETS}{$set}{$gene[$i]}->{TXEND} +
                $self->{INFO}{UPSTREAM} -
                $self->{SETS}{$set}{$gene[$i]}->{NONREPEAT}[$z]->{END};
                
                my $rec_length = $rec_end - $rec_start;
                my $line_height = $position + 2;
                my $line_start = $scale * $rec_start;
                my $line_length = $scale * $rec_length;
                if ($line_start > 500){next};
                if ($line_start + $line_length > 500){$line_length = 500 - $line_start};
                
                $writer->emptyTag(
                        'rect',
                        x => $line_start,
                        y  => $line_height,
                        width => $line_length,
                        height  => '8',
                        rx => $radius,
                        style => 'fill: black;stroke-width: 1; stroke: black;'
                ); 
            };
            
            for my $z (0 .. $#{$self->{SETS}{$set}{$gene[$i]}->{CONSERVED}}){
                
                my $rec_end = $self->{SETS}{$set}{$gene[$i]}->{TXEND} +
                $self->{INFO}{UPSTREAM} -
                $self->{SETS}{$set}{$gene[$i]}->{CONSERVED}[$z]->{START};
                
                my $rec_start =  $self->{SETS}{$set}{$gene[$i]}->{TXEND} +
                $self->{INFO}{UPSTREAM} -
                $self->{SETS}{$set}{$gene[$i]}->{CONSERVED}[$z]->{END};
                
                my $rec_length = $rec_end - $rec_start;
                my $line_height = $position + 2;
                my $line_start = $scale * $rec_start;
                my $line_length = $scale * $rec_length;
                if ($line_start > 500){next};
                if ($line_start + $line_length > 500){$line_length = 500 - $line_start};
                    
                $writer->emptyTag(
                    'rect',
                    x => $line_start,
                    y  => $position,
                    width => $line_length,
                    height  => '12',
                    rx => $radius,
                    style => 'fill: red;stroke-width: 1; stroke: black;'
                ); 
            };
            $position = $position + 30;
        };
        my $ticks = ($self->{INFO}{UPSTREAM} + $self->{INFO}{DOWNSTREAM})/10;
        my $offset = $position + 6;
        $writer->emptyTag(
            'line',
            x1 => '0',
            y1  => $offset,
            x2 => '500',
            y2  => $offset,
            style => 'stroke-width: 2; stroke: black;'
        );
        for my $i (0 .. 10){
            my $tickpos = $i * 50;
            $writer->emptyTag(
                'rect',
                x => $tickpos,
                y  => $offset,
                width => '1',
                height  => '6',
                rx => $radius,
                style => 'fill: black;stroke-width: 1; stroke: black;'
            );
            my $pos = $position + 25;
            my $location = ($i*$ticks) - $self->{INFO}{UPSTREAM};
            $writer->startTag('text',
                x=> $tickpos,
                y=> $pos,
                style => 'font-family:Arial-BoldMT;font-size:10;fill:black;text-anchor:middle;'
            );
            $writer->characters("$location");
            $writer->endTag('text');
        };
        $writer->endTag('g');
        $writer->endTag('svg');
        $writer->end();
        $io->close;
    };
};

sub upstream_factor_view
{
    my ($self, $set, $region) = @_;    
    my $count = 0;
    my $radius = 1;
    my $image_width = 1200;
    my @factors = keys %{$self->{FACTORS}{LIST}};
    my $number = $#factors;
    my $image_height = 50 * $number + 100;

    my $scale = 500/($self->{INFO}{UPSTREAM} + $self->{INFO}{DOWNSTREAM});
    my @gene = keys %{$self->{SEARCHSET}};
    @gene = sort { ($a =~ /(\d+)/)[0] <=> ($b =~ /(\d+)/)[0] } @gene;
    
    for my $i (0 .. $#gene){
        my $file = "$self->{DIRECTORIES}{WORK}\/$set$region$gene[$i]\.svg";
        my $io = new IO::File;
        $io->open ("> $file");
        my $writer = new XML::Writer(OUTPUT => $io, DATA_MODE => 1, DATA_INDENT => 1);
        
        $writer->startTag(
            'svg',
            height => $image_height,
            width  => $image_width
        );
        
        $writer->emptyTag(
            'rect',
            height => $image_height,
            width  => $image_width,
            fill   => 'white'
        );
        
        $writer->startTag(
            'g',
            id => 'motif-page',
            transform => 'translate(100,100)'
        );
    
        my $position = 0;
        if ($self->{SETS}{$set}{$gene[$i]}->{STRAND} =~ /\+/){
            for my $m (0 .. $#factors){
                for my $p (0 .. $#{$self->{SETS}{$set}{$gene[$i]}->{FACTOR}{$factors[$m]}}){
                    my $rec_start = $self->{SETS}{$set}{$gene[$i]}->{FACTOR}{$factors[$m]}[$p]->{START} -
                    $self->{SETS}{$set}{$gene[$i]}->{TXSTART} +
                    $self->{INFO}{UPSTREAM};
                
                    my $rec_end = $self->{SETS}{$set}{$gene[$i]}->{FACTOR}{$factors[$m]}[$p]->{END} -
                    $self->{SETS}{$set}{$gene[$i]}->{TXSTART} +
                    $self->{INFO}{UPSTREAM};
                
                    my $rec_length = $rec_end - $rec_start;
                    my $line_height = $position + 2;
                    my $line_start = $scale * $rec_start;
                    my $line_length = $scale * $rec_length;
                    if ($line_start > 500){next};
                    if ($line_start + $line_length > 500){$line_length = 500 - $line_start};
                
                    $writer->emptyTag(
                        'rect',
                        x => $line_start,
                        y  => $line_height,
                        width => $line_length,
                        height  => '8',
                        style => 'fill: black;stroke-width: 0.1; stroke: black;'
                    );
                };
                my $pos_text = $position + 10;
                $writer->startTag('text',
                    x=> '525',
                    y=> $pos_text,
                    style => 'font-family:Arial-BoldMT;font-size:12;fill:black;text-anchor:front;'
                );
                my $motnum = $m + 1;
                $writer->characters("$factors[$m]");
                $writer->endTag('text');
                $position += 30;
            };
            
            my $offset = $position + 6;
            
            $writer->emptyTag(
                'line',
                x1 => '0',
                y1  => $offset,
                x2 => '500',
                y2  => $offset,
                style => 'stroke-width: 2; stroke: black;'
            );
            
            my $pos_text = $position + 10;
            $writer->startTag('text',
                x=> '525',
                y=> $pos_text,
                style => 'font-family:Arial-BoldMT;font-size:12;fill:black;text-anchor:front;'
            );
            $writer->characters("$gene[$i]");
            $writer->endTag('text');
            
            for my $z (0 .. $#{$self->{SETS}{$set}{$gene[$i]}->{NONREPEAT}}){
                my $rec_start = $self->{SETS}{$set}{$gene[$i]}->{NONREPEAT}[$z]->{START} -
                $self->{SETS}{$set}{$gene[$i]}->{TXSTART} +
                $self->{INFO}{UPSTREAM};
                
                my $rec_end = $self->{SETS}{$set}{$gene[$i]}->{NONREPEAT}[$z]->{END} -
                $self->{SETS}{$set}{$gene[$i]}->{TXSTART} +
                $self->{INFO}{UPSTREAM};
                
                my $rec_length = $rec_end - $rec_start;
                my $line_height = $position + 2;
                my $line_start = $scale * $rec_start;
                my $line_length = $scale * $rec_length;
                if ($line_start > 500){next};
                if ($line_start + $line_length > 500){$line_length = 500 - $line_start};
                
                $writer->emptyTag(
                        'rect',
                        x => $line_start,
                        y  => $line_height,
                        width => $line_length,
                        height  => '8',
                        rx => $radius,
                        style => 'fill: black;stroke-width: 1; stroke: black;'
                ); 
            };
            
            for my $z (0 .. $#{$self->{SETS}{$set}{$gene[$i]}->{CONSERVED}}){
                
                my $rec_start = $self->{SETS}{$set}{$gene[$i]}->{CONSERVED}[$z]->{START} -
                $self->{SETS}{$set}{$gene[$i]}->{TXSTART} +
                $self->{INFO}{UPSTREAM};
                
                my $rec_end = $self->{SETS}{$set}{$gene[$i]}->{CONSERVED}[$z]->{END} -
                $self->{SETS}{$set}{$gene[$i]}->{TXSTART} +
                $self->{INFO}{UPSTREAM};
                
                my $rec_length = $rec_end - $rec_start;
                my $line_height = $position + 2;
                my $line_start = $scale * $rec_start;
                my $line_length = $scale * $rec_length;
                if ($line_start > 500){next};
                if ($line_start + $line_length > 500){$line_length = 500 - $line_start};
                    
                $writer->emptyTag(
                    'rect',
                    x => $line_start,
                    y  => $position,
                    width => $line_length,
                    height  => '12',
                    rx => $radius,
                    style => 'fill: red;stroke-width: 1; stroke: black;'
                ); 
            };
            $position = $position + 30;
        }
        else{
            #add scans here
            for my $m (0 .. $#factors){
                for my $p (0 .. $#{$self->{SETS}{$set}{$gene[$i]}->{FACTORS}{$factors[$m]}}){
                    my $rec_start = $self->{SETS}{$set}{$gene[$i]}->{TXEND} +
                    $self->{INFO}{UPSTREAM} -
                    $self->{SETS}{$set}{$gene[$i]}->{FACTORS}{$factors[$m]}[$p]->{END};
                
                    my $rec_end = 
                    $self->{SETS}{$set}{$gene[$i]}->{TXEND} +
                    $self->{INFO}{UPSTREAM} -
                    $self->{SETS}{$set}{$gene[$i]}->{FACTORS}{$factors[$m]}[$p]->{START};
                
                    my $rec_length = $rec_end - $rec_start;
                    my $line_height = $position + 2;
                    my $line_start = $scale * $rec_start;
                    my $line_length = $scale * $rec_length;
                    if ($line_start > 500){next};
                    if ($line_start + $line_length > 500){$line_length = 500 - $line_start};
                
                    $writer->emptyTag(
                        'rect',
                        x => $line_start,
                        y  => $line_height,
                        width => $line_length,
                        height  => '8',
                        style => 'fill: black;stroke-width: 0.1; stroke: black;'
                    ); 
                };
                my $pos_text = $position + 10;
                $writer->startTag('text',
                    x=> '525',
                    y=> $pos_text,
                    style => 'font-family:Arial-BoldMT;font-size:12;fill:black;text-anchor:front;'
                );
                my $motnum = $m + 1;
                $writer->characters("Motif $motnum");
                $writer->endTag('text');                
                $position += 30;
            };
            
            my $offset = $position + 6;
            
            $writer->emptyTag(
                'line',
                x1 => '0',
                y1  => $offset,
                x2 => '500',
                y2  => $offset,
                style => 'stroke-width: 2; stroke: black;'
            );
            
            my $pos_text = $position + 10;
            $writer->startTag('text',
                x=> '525',
                y=> $pos_text,
                style => 'font-family:Arial-BoldMT;font-size:12;fill:black;text-anchor:front;'
            );
            $writer->characters("$gene[$i]");
            $writer->endTag('text');
            
            for my $z (0 .. $#{$self->{SETS}{$set}{$gene[$i]}->{NONREPEAT}}){
                my $rec_end =  $self->{SETS}{$set}{$gene[$i]}->{TXEND} +
                $self->{INFO}{UPSTREAM} -
                $self->{SETS}{$set}{$gene[$i]}->{NONREPEAT}[$z]->{START};
                
                my $rec_start = $self->{SETS}{$set}{$gene[$i]}->{TXEND} +
                $self->{INFO}{UPSTREAM} -
                $self->{SETS}{$set}{$gene[$i]}->{NONREPEAT}[$z]->{END};
                
                my $rec_length = $rec_end - $rec_start;
                my $line_height = $position + 2;
                my $line_start = $scale * $rec_start;
                my $line_length = $scale * $rec_length;
                if ($line_start > 500){next};
                if ($line_start + $line_length > 500){$line_length = 500 - $line_start};
                
                $writer->emptyTag(
                        'rect',
                        x => $line_start,
                        y  => $line_height,
                        width => $line_length,
                        height  => '8',
                        rx => $radius,
                        style => 'fill: black;stroke-width: 1; stroke: black;'
                ); 
            };
            
            for my $z (0 .. $#{$self->{SETS}{$set}{$gene[$i]}->{CONSERVED}}){
                
                my $rec_end = $self->{SETS}{$set}{$gene[$i]}->{TXEND} +
                $self->{INFO}{UPSTREAM} -
                $self->{SETS}{$set}{$gene[$i]}->{CONSERVED}[$z]->{START};
                
                my $rec_start =  $self->{SETS}{$set}{$gene[$i]}->{TXEND} +
                $self->{INFO}{UPSTREAM} -
                $self->{SETS}{$set}{$gene[$i]}->{CONSERVED}[$z]->{END};
                
                my $rec_length = $rec_end - $rec_start;
                my $line_height = $position + 2;
                my $line_start = $scale * $rec_start;
                my $line_length = $scale * $rec_length;
                if ($line_start > 500){next};
                if ($line_start + $line_length > 500){$line_length = 500 - $line_start};
                    
                $writer->emptyTag(
                    'rect',
                    x => $line_start,
                    y  => $position,
                    width => $line_length,
                    height  => '12',
                    rx => $radius,
                    style => 'fill: red;stroke-width: 1; stroke: black;'
                ); 
            };
            $position = $position + 30;
        };
        my $ticks = ($self->{INFO}{UPSTREAM} + $self->{INFO}{DOWNSTREAM})/10;
        my $offset = $position + 6;
        $writer->emptyTag(
            'line',
            x1 => '0',
            y1  => $offset,
            x2 => '500',
            y2  => $offset,
            style => 'stroke-width: 2; stroke: black;'
        );
        for my $i (0 .. 10){
            my $tickpos = $i * 50;
            $writer->emptyTag(
                'rect',
                x => $tickpos,
                y  => $offset,
                width => '1',
                height  => '6',
                rx => $radius,
                style => 'fill: black;stroke-width: 1; stroke: black;'
            );
            my $pos = $position + 25;
            my $location = ($i*$ticks) - $self->{INFO}{UPSTREAM};
            $writer->startTag('text',
                x=> $tickpos,
                y=> $pos,
                style => 'font-family:Arial-BoldMT;font-size:10;fill:black;text-anchor:middle;'
            );
            $writer->characters("$location");
            $writer->endTag('text');
        };
        $writer->endTag('g');
        $writer->endTag('svg');
        $writer->end();
        $io->close;
    };
};

sub upstream_sequence_view
{
    my ($self, $set, $file) = @_;    
    my $count = 0;
    my $radius = 1;
    
    my $image_width = 1200;
    my $scale = 500/($self->{INFO}{UPSTREAM} + $self->{INFO}{DOWNSTREAM});
    my @gene = keys %{$self->{SETS}{$set}};
    my $image_height = 50 * $#gene;
    
    @gene = sort { ($a =~ /(\d+)/)[0] <=> ($b =~ /(\d+)/)[0] } @gene;
    my $io = new IO::File;
    $io->open ("> $file");
    my $writer = new XML::Writer(OUTPUT => $io, DATA_MODE => 1, DATA_INDENT => 1);
    $writer->startTag(
        'svg',
        height => $image_height,
        width  => $image_width
    );
        
    $writer->emptyTag(
        'rect',
        height => $image_height,
        width  => $image_width,
        fill   => 'white'
    );
        
    $writer->startTag(
        'g',
        id => 'motif-page',
        transform => 'translate(100,100)'
    );
    my $position = 0;
    for my $i (0 .. $#gene){
        if ($self->{SETS}{$set}{$gene[$i]}->{STRAND} =~ /\+/)
        {
            my $offset = $position + 6;
            $writer->emptyTag(
                'line',
                x1 => '0',
                y1  => $offset,
                x2 => '500',
                y2  => $offset,
                style => 'stroke-width: 2; stroke: black;'
            );
            
            my $pos_text = $position + 10;
            $writer->startTag('text',
                x=> '525',
                y=> $pos_text,
                style => 'font-family:Arial-BoldMT;font-size:12;fill:black;text-anchor:front;'
            );
            $writer->characters("$gene[$i]");
            $writer->endTag('text');
            
            for my $z (0 .. $#{$self->{SETS}{$set}{$gene[$i]}->{NONREPEAT}}){
                my $rec_start = $self->{SETS}{$set}{$gene[$i]}->{NONREPEAT}[$z]->{START} -
                $self->{SETS}{$set}{$gene[$i]}->{TXSTART} +
                $self->{INFO}{UPSTREAM};
                
                my $rec_end = $self->{SETS}{$set}{$gene[$i]}->{NONREPEAT}[$z]->{END} -
                $self->{SETS}{$set}{$gene[$i]}->{TXSTART} +
                $self->{INFO}{UPSTREAM};
                
                my $rec_length = $rec_end - $rec_start;
                my $line_height = $position + 2;
                my $line_start = $scale * $rec_start;
                my $line_length = $scale * $rec_length;
                if ($line_start > 500){next};
                if ($line_start + $line_length > 500){$line_length = 500 - $line_start};
                
                $writer->emptyTag(
                        'rect',
                        x => $line_start,
                        y  => $line_height,
                        width => $line_length,
                        height  => '8',
                        rx => $radius,
                        style => 'fill: black;stroke-width: 1; stroke: black;'
                ); 
            };
            
            for my $z (0 .. $#{$self->{SETS}{$set}{$gene[$i]}->{CONSERVED}}){
                
                my $rec_start = $self->{SETS}{$set}{$gene[$i]}->{CONSERVED}[$z]->{START} -
                $self->{SETS}{$set}{$gene[$i]}->{TXSTART} +
                $self->{INFO}{UPSTREAM};
                
                my $rec_end = $self->{SETS}{$set}{$gene[$i]}->{CONSERVED}[$z]->{END} -
                $self->{SETS}{$set}{$gene[$i]}->{TXSTART} +
                $self->{INFO}{UPSTREAM};
                
                my $rec_length = $rec_end - $rec_start;
                my $line_height = $position + 2;
                my $line_start = $scale * $rec_start;
                my $line_length = $scale * $rec_length;
                if ($line_start > 500){next};
                if ($line_start + $line_length > 500){$line_length = 500 - $line_start};
                    
                $writer->emptyTag(
                    'rect',
                    x => $line_start,
                    y  => $position,
                    width => $line_length,
                    height  => '12',
                    rx => $radius,
                    style => 'fill: red;stroke-width: 1; stroke: black;'
                ); 
            };
            $position = $position + 30;
        }
        else{
            my $offset = $position + 6;
            $writer->emptyTag(
                'line',
                x1 => '0',
                y1  => $offset,
                x2 => '500',
                y2  => $offset,
                style => 'stroke-width: 2; stroke: black;'
            );
            
            my $pos_text = $position + 10;
            $writer->startTag('text',
                x=> '525',
                y=> $pos_text,
                style => 'font-family:Arial-BoldMT;font-size:12;fill:black;text-anchor:front;'
            );
            $writer->characters("$gene[$i]");
            $writer->endTag('text');
            
            for my $z (0 .. $#{$self->{SETS}{$set}{$gene[$i]}->{NONREPEAT}}){
                my $rec_end =  $self->{SETS}{$set}{$gene[$i]}->{TXEND} +
                $self->{INFO}{UPSTREAM} -
                $self->{SETS}{$set}{$gene[$i]}->{NONREPEAT}[$z]->{START};
                
                my $rec_start = $self->{SETS}{$set}{$gene[$i]}->{TXEND} +
                $self->{INFO}{UPSTREAM} -
                $self->{SETS}{$set}{$gene[$i]}->{NONREPEAT}[$z]->{END};
                
                my $rec_length = $rec_end - $rec_start;
                my $line_height = $position + 2;
                my $line_start = $scale * $rec_start;
                my $line_length = $scale * $rec_length;
                if ($line_start > 500){next};
                if ($line_start + $line_length > 500){$line_length = 500 - $line_start};
                
                $writer->emptyTag(
                        'rect',
                        x => $line_start,
                        y  => $line_height,
                        width => $line_length,
                        height  => '8',
                        rx => $radius,
                        style => 'fill: black;stroke-width: 1; stroke: black;'
                ); 
            };
            
            for my $z (0 .. $#{$self->{SETS}{$set}{$gene[$i]}->{CONSERVED}}){
                
                my $rec_end = $self->{SETS}{$set}{$gene[$i]}->{TXEND} +
                $self->{INFO}{UPSTREAM} -
                $self->{SETS}{$set}{$gene[$i]}->{CONSERVED}[$z]->{START};
                
                my $rec_start =  $self->{SETS}{$set}{$gene[$i]}->{TXEND} +
                $self->{INFO}{UPSTREAM} -
                $self->{SETS}{$set}{$gene[$i]}->{CONSERVED}[$z]->{END};
                
                my $rec_length = $rec_end - $rec_start;
                my $line_height = $position + 2;
                my $line_start = $scale * $rec_start;
                my $line_length = $scale * $rec_length;
                if ($line_start > 500){next};
                if ($line_start + $line_length > 500){$line_length = 500 - $line_start};
                    
                $writer->emptyTag(
                    'rect',
                    x => $line_start,
                    y  => $position,
                    width => $line_length,
                    height  => '12',
                    rx => $radius,
                    style => 'fill: red;stroke-width: 1; stroke: black;'
                ); 
            };
            $position = $position + 30;
        };
    };
    my $ticks = ($self->{INFO}{UPSTREAM} + $self->{INFO}{DOWNSTREAM})/10;
    my $offset = $position + 6;
    $writer->emptyTag(
        'line',
        x1 => '0',
        y1  => $offset,
        x2 => '500',
        y2  => $offset,
        style => 'stroke-width: 2; stroke: black;'
    );
    for my $i (0 .. 10){
        my $tickpos = $i * 50;
        $writer->emptyTag(
            'rect',
            x => $tickpos,
            y  => $offset,
            width => '1',
            height  => '6',
            rx => $radius,
            style => 'fill: black;stroke-width: 1; stroke: black;'
        );
        my $pos = $position + 25;
        my $location = ($i*$ticks) - $self->{INFO}{UPSTREAM};
        $writer->startTag('text',
            x=> $tickpos,
            y=> $pos,
            style => 'font-family:Arial-BoldMT;font-size:10;fill:black;text-anchor:middle;'
        );
        $writer->characters("$location");
        $writer->endTag('text');
    };
    $writer->endTag('g');
    $writer->endTag('svg');
    $writer->end();
    $io->close;
};

1;
