#!/usr/bin/perl -w
#
# find_anomalies_public_config.pl
#
# by Gary Williams
#
# Aug 2008
#
###################################################################################
# This is the general, public version of the system used in the WormBase project to
# find curation anomalies run from a config file.
# It has had many worm-spcific features removed.
###################################################################################
#
# This looks for anomalous things such as protein homologies not
# matching a CDS and stores the results in in a data file ready to be read into the SQL database
#
# Last updated by: $Author: gw3 $
# Last updated on: $Date: 2008/08/26 12:28:48 $

# Problems and assumptions with making this generalised:
#
# We assume that data is available as GFF files
#
# We assume that the genomic sequence is available as FASTA-format
#files
#
# We assume that the BLAST score is available in the GFF
#
# We assume that the region of homology is specified in the GFF by the
# start and end position after the ID name in the last field of the
# GFF. This is the standard in the GFF3 'Target=ID start end' part of
# the attribute field.
#
# We don't know the orientation of the EST/mRNA/etc reads. (all
# comparisons with these have to be made sense-independent)
#
# We don't know how isoforms are to be distinguished from genes coding
# for a single protein. (Have a routine that is re-written for
# specific systems)
#
# We don't know how to convert isoform IDs to a standard CDS ID (Have
# a routine that is re-written for specific systems)
#
# We don't know how to convert gene transcript IDs to a standard CDS
# ID (Have a routine that is re-written for specific systems)
#
# We don't know what types of EST/mRNA/RST/OST/etc. transcripts are
# available (will have to have a configuration file)
#
# We don't know what types of gene predictions are available that
# differ from the gene models. (will have to change the program somewhat)
#
# We have used a Position Weight Matrix derived from C. elegans splice
# sites. The new target species will have to haev one set up for it.
# (Need to count frequencies of nucleotides in the region of known
# splice sites)
#
# For each chromosome, the GFF data and the chromosamal sequence are
# all held in memory - this is trivial to do with nematodes, but may
# be a problem for larger genomes. (May need a rewrite)

use strict;

#use lib $ENV{'CVS_DIR'};
use Getopt::Long;
use Carp;
use Overlap_config;
use PWM;
use POSIX qw(log10);

######################################
# variables and command-line options #
######################################

my ( $help, $verbose, $config, $output );

GetOptions(
    "help"     => \$help,
    "verbose"  => \$verbose,
    "config:s" => \$config,    # file to specify lots of things for your system
    "output:s" => \$output,    # output file to be put into SQL database
);

# Display help if required
&usage("Help") if ($help);

# check on the output file
if ( !defined $output ) { die "-output is not specified\n"; }

# open the config file and read
# species, database, overlap config file
# lists of GFF data, list of chromosomes
# list of which analyses are being done by which species
if ( !defined $config ) {
    $config = "./find_anomalies.config";
    print "Using config file: $config\n";
}
my %config = &read_config($config);

# get some useful data from the config hash
my $species = $config{species};    # name of the species
if ( scalar @{ $config{run} } == 0 ) {
    die "There are no analyses specified to be RUN in the config file.\n";
}
my @run            = @{ $config{run} };            # list of which analyses to run
my @chromosomes    = @{ $config{chromosomes} };    # list of chromosomes
my $overlap_config = $config{overlap};        # config file for the Overlap module
my $dnadir         = $config{dnadir};         # directory containing the chromosome FASTA sequence file(s)
my %define         = %{ $config{define} };    # set of definitions of the analyses

# global variable holding the GFF data
my %GFF_data;

##########################
# MAIN BODY OF SCRIPT
##########################

my %dna_entry;   # store for dna sequences for when reading multi-entry fasta sequence file

my $ovlp;    # Overlap object

# output file of data to write to database
open( DAT, "> $output" ) || die "Can't open $output: $!\n";

# and output the species line for the data file
print DAT "SPECIES\t$species\n\n";

my $pwm = PWM->new;

# hash of anomaly counts
my %anomaly_count;

# now write the data commands to delete all the previous analyses from
# the SQL database
# go through the definitions of the analyses to run finding the names of the anomalies
# these are in the defnition like 'output <GFF_LIST> as <ANOMALY_NAME>'
foreach my $analysis (@run) {
    my $def = $define{$analysis};
    if ( !defined $def ) {
        die "ERROR: The analysis $analysis has not been defined\n";
    }

    # get each of the anomaly names output from the definition of this analysis
    while ( $def =~ /OUTPUT\s+\S+\s+AS\s+(\S+)/ig ) {
        &delete_anomalies($1);
    }
}

######################################################################
foreach my $chromosome (@chromosomes) {
    print "Processing chromosome $chromosome\n" if ($verbose);

    # delete any GFF lists from previous chromosomes
    %GFF_data = ();

    # get the Overlap object
    $ovlp = Overlap_config->new( glob($overlap_config) );

    print "Finding anomalies\n\n" if ($verbose);

    foreach my $analysis (@run) {
        my $def = $define{$analysis};
        print "Running analysis $analysis\n" if ($verbose);
        &parse_analysis( $chromosome, $analysis, $def );
    }
}
######################################################################

# close the output datafile
close(DAT);

# report the numbers of anomalies found
print "\nCount of anomalies found:\n";
foreach my $anomaly_type ( keys %anomaly_count ) {
    my $count = $anomaly_count{$anomaly_type};
    print "$anomaly_type\t$count\n";
}

print "Finished.\n" if ($verbose);
exit(0);

##############################################################
#
# Subroutines
#
##############################################################

##########################################
#
# Read the config file
# my %config = &read_config($config);

sub read_config {
    my $config_file = shift;
    my %config;

    open( CONFIG, "<$config_file" ) || die "Can't open file $config_file: $!\n";

    while ( my $line = <CONFIG> ) {
        chomp $line;

        if ( $line =~ /^\s*$/ || $line =~ /^\s*#/ ) {    # skip blank lines and comments
            next;
        }
        elsif ( $line =~ /^species\s+(\S+)/i ) {
            # name of the species
            $config{species} = $1;
        }
        elsif ( $line =~ /^run\s+(\S+)/i ) {
            # one or more lines of which analyses to run
            push @{ $config{run} }, $1;
        }
        elsif ( $line =~ /^chromosomes\s+(.+)/i ) {
            # one or more lines of lists of chromosomes separated by commas or spaces
            push @{ $config{chromosomes} }, split /\s+/, $1;
        }
        elsif ( $line =~ /^overlap\s+(\S+)/i ) {
            # config file for the Overlap module
            $config{overlap} = $1;
        }
        elsif ( $line =~ /^dnadir\s+(\S+)/i ) {
            # the directory containing the chromosomal sequence files
            $config{dnadir} = $1;
        }
        elsif ( $line =~ /^define\s+(\S+)\s+\{/i ) {
            while ( my $line2 = <CONFIG> ) {
                chomp $line2;

                if ( $line2 =~ /^\s*$/ || $line2 =~ /^\s*#/ ) {
                    next;
                }    # skip blank lines and comments

                if ( $line2 =~ /^\s*\}\s*$/ ) {
                    last;
                }    # a '}' on its own ends the definition
                $config{define}{$1} .= $line2 . " ";
            }
        }
        else {
            die "Unrecognised line in $config_file:\n$line\n";
        }
    }
    close(CONFIG);

    return %config;
}

##########################################
# parse_analysis
# takes a definition of the analysis from the config file and runs the
# analysis based on the defniition

sub parse_analysis {
    my ( $chromosome, $analysis, $def ) = @_;

    # the state of this analysis
    my %state;

    my @cmds = split /(MINSCORE\s+|MATCH\s+|OUTPUT\s+|MERGE\s+|FIRST_HSP\s+|TRUNCATE\s+|DROP\s+|ID_PATTERN\s+|MAXLENGTH\s+|MINLENGTH\s+|ISOFORM\s+|COUNT_ID\s+|READ_GFF\s+|CONFIRMED_INTRONS\s+|SPLIT_ENDS\s+|SPLICE_SCORE\s+|FRAMESHIFT\s+|PROTEIN_ALIGN\s+)/i, $def;

# e.g.
# 0  '  '
# 1  'MATCH '
# 2  'GFF_name1 {     VS GFF_name2 OTHER_SENSE   } => match_name, mismatch_name   '
# 3  'OUTPUT '
# 4  'match_name AS anomaly_name SCORE 1 '

    for ( my $i = 0 ; $i < scalar @cmds ; $i++ ) {
        if ( $cmds[$i] =~ /^\s*$/ ) {
            next;
        }
        elsif ( $cmds[$i] =~ /OUTPUT/i ) {
            &output( $chromosome, $analysis, $cmds[ ++$i ], \%state );
        }
        elsif ( $cmds[$i] =~ /MINSCORE/i ) {
            &minscore( $chromosome, $analysis, $cmds[ ++$i ], \%state );
        }
        elsif ( $cmds[$i] =~ /MINLENGTH/i ) {
            &minlength( $chromosome, $analysis, $cmds[ ++$i ], \%state );
        }
        elsif ( $cmds[$i] =~ /MAXLENGTH/i ) {
            &maxlength( $chromosome, $analysis, $cmds[ ++$i ], \%state );
        }
        elsif ( $cmds[$i] =~ /MATCH/i ) {
            &match( $chromosome, $analysis, $cmds[ ++$i ], \%state );
        }
        elsif ( $cmds[$i] =~ /MERGE/i ) {
            &merge( $chromosome, $analysis, $cmds[ ++$i ], \%state );
        }
        elsif ( $cmds[$i] =~ /FIRST_HSP/i ) {
            &first_hsp( $chromosome, $analysis, $cmds[ ++$i ], \%state );
        }
        elsif ( $cmds[$i] =~ /TRUNCATE/i ) {
            &truncate( $chromosome, $analysis, $cmds[ ++$i ], \%state );
        }
        elsif ( $cmds[$i] =~ /DROP/i ) {
            &drop( $chromosome, $analysis, $cmds[ ++$i ], \%state );
        }
        elsif ( $cmds[$i] =~ /ID_PATTERN/i ) {
            &id_pattern( $chromosome, $analysis, $cmds[ ++$i ], \%state );
        }
        elsif ( $cmds[$i] =~ /ISOFORM/i ) {
            &isoform( $chromosome, $analysis, $cmds[ ++$i ], \%state );
        }
        elsif ( $cmds[$i] =~ /COUNT_ID/i ) {
            &count_id( $chromosome, $analysis, $cmds[ ++$i ], \%state );
        }
        elsif ( $cmds[$i] =~ /READ_GFF/i ) {
            &read_gff( $chromosome, $analysis, $cmds[ ++$i ], \%state );
        }
        elsif ( $cmds[$i] =~ /CONFIRMED_INTRONS/i ) {
            &confirmed_introns( $chromosome, $analysis, $cmds[ ++$i ], \%state );
        }
        elsif ( $cmds[$i] =~ /SPLIT_ENDS/i ) {
            &split_ends( $chromosome, $analysis, $cmds[ ++$i ], \%state );
        }
        elsif ( $cmds[$i] =~ /SPLICE_SCORE/i ) {
            &splice_score( $chromosome, $analysis, $cmds[ ++$i ], \%state );
        }
        elsif ( $cmds[$i] =~ /FRAMESHIFT/i ) {
            &frameshift( $chromosome, $analysis, $cmds[ ++$i ], \%state );
        }
        elsif ( $cmds[$i] =~ /PROTEIN_ALIGN/i ) {
            &protein_align( $chromosome, $analysis, $cmds[ ++$i ], \%state );
        }
        else {
            die "ERROR: Don't understand the command '$cmds[$i] $cmds[$i+1]' in the definition of analysis $analysis\n";
        }
    }
}

##########################################
# remove features less than a minimum score
sub minscore {
    my ( $chromosome, $analysis, $cmd, $state_href ) = @_;

    print "chr=$chromosome analysis=$analysis cmd=MINSCORE $cmd\n" if ($verbose);

    my ( $gff1_name, $threshold, $gff2_name ) = ( $cmd =~ /(\S+)\s+([+-]*\d+)\s*=>\s*(\S+)/i );
    my $gff1_array = &get_gff( $chromosome, $gff1_name );
    my @new_gff;

    foreach my $feature ( @{$gff1_array} ) {    # id, $chrom_start, $chrom_end, $chrom_strand
        my $score = $feature->[6];
        if ( $score >= $threshold ) {
            push @new_gff, $feature;
        }
    }

    @new_gff = sort { $a->[1] <=> $b->[1] or $a->[2] <=> $b->[2] } @new_gff;

    # store it
    $GFF_data{$gff2_name} = \@new_gff;

    print "Have ", scalar @new_gff, " MINSCORE\n" if ($verbose);

}

##########################################
# remove features less than a minimum length
sub minlength {
    my ( $chromosome, $analysis, $cmd, $state_href ) = @_;

    print "chr=$chromosome analysis=$analysis cmd=MINLENGTH $cmd\n" if ($verbose);

    my ( $gff1_name, $threshold, $gff2_name ) = ( $cmd =~ /(\S+)\s+([+-]*\d+)\s*=>\s*(\S+)/i );
    my $gff1_array = &get_gff( $chromosome, $gff1_name );
    my @new_gff;

    foreach my $feature ( @{$gff1_array} ) {    # id, $chrom_start, $chrom_end, $chrom_strand
        my $start = $feature->[1];
        my $end   = $feature->[2];
        if ( $end - $start + 1 >= $threshold ) {
            push @new_gff, $feature;
        }
    }

    @new_gff = sort { $a->[1] <=> $b->[1] or $a->[2] <=> $b->[2] } @new_gff;

    # store it
    $GFF_data{$gff2_name} = \@new_gff;

    print "Have ", scalar @new_gff, " MINLENGTH\n" if ($verbose);

}

##########################################
# remove features more than a maximum length
sub maxlength {
    my ( $chromosome, $analysis, $cmd, $state_href ) = @_;

    print "chr=$chromosome analysis=$analysis cmd=MAXLENGTH $cmd\n" if ($verbose);

    my ( $gff1_name, $threshold, $gff2_name ) = ( $cmd =~ /(\S+)\s+([+-]*\d+)\s*=>\s*(\S+)/i );
    my $gff1_array = &get_gff( $chromosome, $gff1_name );
    my @new_gff;

    foreach my $feature ( @{$gff1_array} ) {    # id, $chrom_start, $chrom_end, $chrom_strand
        my $start = $feature->[1];
        my $end   = $feature->[2];
        if ( $end - $start + 1 <= $threshold ) {
            push @new_gff, $feature;
        }
    }

    @new_gff = sort { $a->[1] <=> $b->[1] or $a->[2] <=> $b->[2] } @new_gff;

    # store it
    $GFF_data{$gff2_name} = \@new_gff;

    print "Have ", scalar @new_gff, " MAXLENGTH\n" if ($verbose);

}

##########################################
# match one or more GFF lists against a main GFF list
sub match {
    my ( $chromosome, $analysis, $cmd, $state_href ) = @_;

    print "chr=$chromosome analysis=$analysis cmd=MATCH $cmd\n" if ($verbose);

    my ( $gff1_name, $other_gffs, $match_result, $mismatch_result ) = ( $cmd =~ /(\S+)\s+\{(.+)\}\s*=>\s*(\S+?)[,\s]+(\S+)/i );

    #print "$gff1_name, $other_gffs, $match_result, $mismatch_result\n";

    # make sure we have the main GFF list
    my $gff1_array = &get_gff( $chromosome, $gff1_name );

    # get the secondary GFF lists
    my @match_objects;
    my @secondary = split /VS\s+/i, $other_gffs;
    my @required_number_of_matches;    # the minimum number of secondary features we must match
    my @complete; # if set for a secondary list, we want the main feature to be completely covered by the secondary feature

    for ( my $i = 0 ; $i < scalar @secondary ; $i++ ) {

        #print "Get details of secondary match: $secondary[$i]\n";
        if ( $secondary[$i] =~ /^\s*$/ ) {
            next;
        }
        else {
            my ( $gff2_name, $args ) = ( $secondary[$i] =~ /(\S+)\s*(.+)/ );

            #print "gff2_name=$gff2_name, args=$args\n";
            # make sure we have the secondary GFF list
            my $gff2_array = &get_gff( $chromosome, $gff2_name );

            my $required_number_of_matches =
              1;    # the minimum number of secondary features we must match
            my $complete = 0
              ; # we don't require the main feature to be completely covered by the secondary feature by default

            # get the match arguments
            my %args = (
                near_5 =>
                  0,    # 0= we don't allow near 5' matches to count as a match
                near_3     => 0,    # ditto for 3'
                same_sense => 0
                , # 0= we want to allow a match to an object in either sense, set to true if want the same sense only
                other_sense => 0
                , # 0= we want to allow a match to an object in either sense, set to true if want the opposite sense only
                exact_match => 0
                , # 0= we do not want an exact match, but some sort of sloppy overlap
            );
            my @args = split /\s+/, $args;

            for ( my $j = 0 ; $j < scalar @args ; $j++ ) {
                if ( $args[$j] =~ /^\s*$/ ) {
                    next;
                }
                elsif ( $args[$j] =~ /NEAR/i ) {
                    $args{near_5} = $args[ ++$j ];
                    $args{near_3} = $args{near_5};
                }
                elsif ( $args[$j] =~ /NEAR_5/i ) {
                    $args{near_5} = $args[ ++$j ];
                }
                elsif ( $args[$j] =~ /NEAR_3/i ) {
                    $args{near_3} = $args[ ++$j ];
                }
                elsif ( $args[$j] =~ /SAME_SENSE/i ) {
                    $args{same_sense} = 1;
                    if ( $args{other_sense} ) {
                        die "ERROR: Can't have SAME_SENSE and OTHER_SENSE in the command '$secondary[$i]' in the definition of analysis $analysis\n";
                    }
                }
                elsif ( $args[$j] =~ /OTHER_SENSE/i ) {
                    $args{other_sense} = 1;
                    if ( $args{same_sense} ) {
                        die "ERROR: Can't have SAME_SENSE and OTHER_SENSE in the command '$secondary[$i]' in the definition of analysis $analysis\n";
                    }
                }
                elsif ( $args[$j] =~ /EXACT/i ) {
                    $args{exact_match} = 1;
                }
                elsif ( $args[$j] =~ /MINIMUM_OVERLAPS/i )
                {   # set the minimum number of secondary features we must match
                    $required_number_of_matches = $args[ ++$j ];
                }
                elsif ( $args[$j] =~ /COMPLETE/i )
                { # we want the main feature to be completely covered by the secondary feature
                    $complete = 1;
                }
                else {
                    die "ERROR: Don't understand '$args[$j]' in the command '$secondary[$i]' in the definition of analysis $analysis\n";
                }
            }

            # make the secondary match object and associated arrays
            my $match_obj = $ovlp->compare( $gff2_array, %args );
            push @match_objects,              $match_obj;
            push @required_number_of_matches, $required_number_of_matches;
            push @complete,                   $complete;
        }
    }

    #
    # run the match
    #
    my ( @matched_gff, @unmatched_gff );

    foreach my $feature ( @{$gff1_array} ) {    # $id, $chrom_start, $chrom_end, $chrom_strand

        my $id          = $feature->[0];
        my $got_a_match = 0;               # not yet seen a match to anything

        my $i = 0;
      SECONDARY_LIST:
        foreach my $match_object (@match_objects) {
            my $required_number_of_matches = $required_number_of_matches[$i];
            my $complete                   = $complete[$i];
            my @matching_gff2              = $match_object->match($feature);
            if ( scalar @matching_gff2 >= $required_number_of_matches ) {

                #print "Got a match for $id\n";

                if ($complete)
                { # see if we require the main feature to be completely covered by the secondary feature
                    foreach my $result (@matching_gff2) {
                        my ( $prop1, $prop2 ) =
                          $match_object->matching_proportions($result);
                        if ( $prop1 == 1 ) {    # complete overlap found
                            $got_a_match = 1;
                            last SECONDARY_LIST
                              ; # don't bother to match against the other secondary GFF lists
                        }
                    }
                }
                else {
                    $got_a_match = 1;
                    last SECONDARY_LIST
                      ; # don't bother to match against the other secondary GFF lists
                }

            }
            $i++;
        }

        # did we have a match or not?
        if ($got_a_match) {
            push @matched_gff, $feature;
        }
        else {
            push @unmatched_gff, $feature;
        }
    }

    # store the results
    $GFF_data{$match_result}    = \@matched_gff;
    $GFF_data{$mismatch_result} = \@unmatched_gff;

    print "Have ", scalar @{ $GFF_data{$match_result} }, " MATCHES and ",
      scalar @{ $GFF_data{$mismatch_result} }, " MISMATCHES\n"
      if ($verbose);

}

##########################################
# output a GFF list as anomalies
sub output {
    my ( $chromosome, $analysis, $cmd, $state_href ) = @_;

    print "chr=$chromosome analysis=$analysis cmd=OUTPUT $cmd\n" if ($verbose);

    my $score = 1;
    my $minscore;
    my $maxscore;
    my $logscore = 0;

    my ( $gff_name, $anomaly, $other ) = ( $cmd =~ /(\S+)\s+AS\s+(\S+)\s*(.*)/i );

    #print "$gff_name, $anomaly, $other\n";

    my @cmds = split /(SCORE\s+|SCOREMIN\s+|SCOREMAX\s+|LOGSCORE\s+)/i, $other;
    for ( my $i = 0 ; $i < scalar @cmds ; $i++ ) {
        if ( $cmds[$i] =~ /^\s*$/ ) {
            next;

        }
        elsif ( $cmds[$i] =~ /SCORE/i ) {    # use this value as the score
            $score = $cmds[ ++$i ];
        }
        elsif ( $cmds[$i] =~ /SCOREMIN/i )
        {    # if score of feature is less than this, set score to this
            $minscore = $cmds[ ++$i ];
        }
        elsif ( $cmds[$i] =~ /SCOREMAX/i )
        {    # if score of feature is more than this, set score to this
            $maxscore = $cmds[ ++$i ];
        }
        elsif ( $cmds[$i] =~ /LOGSCORE/i )
        {    # convert score of feature to LOG10 and output that value
            $logscore = 1;
        }
        else {
            die "ERROR: Don't understand the command '$cmds[$i]' in the definition of analysis $analysis\n";
        }
    }

    # the GFF_list should exist by now
    if ( !exists $GFF_data{$gff_name} ) {

        # list the existing GFF lists
        print "ERROR: In analysis '$analysis' the result GFF list '$gff_name' does not exist\n";
        print "Existing GFF lists:\n";
        foreach my $key ( keys %GFF_data ) {
            print "'$key'\n";
        }
        die "\n";
    }

    foreach my $feature ( @{ $GFF_data{$gff_name} } ) {
        my $id            = $feature->[0];
        my $start         = $feature->[1];
        my $end           = $feature->[2];
        my $strand        = $feature->[3];
        my $explanation   = '';             # to be done +++
        my $anomaly_score = $feature->[6];  # use the existing score, if present

        if ( defined $anomaly_score && defined $logscore ) {
            if ( $anomaly_score < 0 ) { $anomaly_score = 0 }
            $anomaly_score = POSIX::log10($anomaly_score);
        }

        if (   defined $anomaly_score
            && defined $minscore
            && $minscore > $anomaly_score )
        {
            $anomaly_score = $minscore;
        }

        if (   defined $anomaly_score
            && defined $maxscore
            && $maxscore < $anomaly_score )
        {
            $anomaly_score = $maxscore;
        }

        if ( defined $score ) { $anomaly_score = $score }

        &output_to_database( $anomaly, $chromosome, $id, $start, $end, $strand, $anomaly_score, $explanation );
    }
    print "Have ", scalar @{ $GFF_data{$gff_name} }, " OUTPUT\n" if ($verbose);

}

##########################################
# merge two or more lists
sub merge {
    my ( $chromosome, $analysis, $cmd, $state_href ) = @_;

    print "chr=$chromosome analysis=$analysis cmd=MERGE $cmd\n" if ($verbose);

    my ( $gff1_name, $operator, $gff2_name, $gff3_name ) =
      ( $cmd =~ /(\S+)\s+(\S+)\s+(\S+)\s*=>\s*(\S+)/i );
    my $gff1_array = &get_gff( $chromosome, $gff1_name );
    my $gff2_array = &get_gff( $chromosome, $gff2_name );

    # operators AND , OR
    # OR  - simple concatenation then sort with no duplicates
    # AND - simple concatenation then sort and only get exact position matches
    my @new_gff;
    my $prev_start  = -1;
    my $prev_end    = -1;
    my $prev_strand = "";
    foreach my $feature (
        sort {
                 $a->[1] <=> $b->[1]
              or $a->[2] <=> $b->[2]
              or $a->[3] cmp $b->[3]
        } ( @{$gff1_array}, @{$gff2_array} )
      )
    {    # id, $chrom_start, $chrom_end, $chrom_strand
        my $start  = $feature->[1];
        my $end    = $feature->[2];
        my $strand = $feature->[3];
        my $same   = 0;
        if (   $start == $prev_start
            && $end == $prev_end
            && $strand eq $prev_strand )
        {
            $same = 1;
        }
        if ( $operator =~ /^OR$/i ) {
            if ( !$same ) {
                push @new_gff, $feature;
            }
        }
        elsif ( $operator =~ /^AND$/i ) {
            if ($same) {
                push @new_gff, $feature;
            }
        }
        else {
            die "ERROR: Don't understand the MERGE operator '$operator' in the definition of analysis $analysis\n";
        }
        $prev_start  = $start;
        $prev_end    = $end;
        $prev_strand = $strand;
    }

    # change the hash to an array sorted by start position
    @new_gff = sort { $a->[1] <=> $b->[1] or $a->[2] <=> $b->[2] } @new_gff;

    # store it
    $GFF_data{$gff3_name} = \@new_gff;

    print "Have ", scalar @new_gff, " MERGE\n" if ($verbose);

}

##########################################
# get just the first (5'-most) HSP of a set of hits for each ID
sub first_hsp {
    my ( $chromosome, $analysis, $cmd, $state_href ) = @_;

    print "chr=$chromosome analysis=$analysis cmd=FIRST_HSP $cmd\n" if ($verbose);

    my ( $gff1_name, $gff2_name ) = ( $cmd =~ /(\S+)\s*=>\s*(\S+)/i );
    my $gff1_array = &get_gff( $chromosome, $gff1_name );

    my %seen_this;

    foreach my $feature ( @{$gff1_array} ) {    # id, $chrom_start, $chrom_end, $chrom_strand
        my $id = $feature->[0];
        if ( $feature->[3] eq '+' && !exists $seen_this{$id} ) {
            $seen_this{$id} = $feature;
        }
        if ( $feature->[3] eq '-' ) {
            $seen_this{$id} = $feature;  # get the last one of these in the list
        }
    }

    # change the hash to an array sorted by start position
    my @new_gff;
    foreach my $id ( keys %seen_this ) {
        push @new_gff, $seen_this{$id};
    }
    @new_gff = sort { $a->[1] <=> $b->[1] or $a->[2] <=> $b->[2] } @new_gff;

    # store it
    $GFF_data{$gff2_name} = \@new_gff;

    print "Have ", scalar @new_gff, " FIRST_HSPs\n" if ($verbose);

}

##########################################
# restrict the features to a region of a few bases before and after the start of the feature
sub truncate {
    my ( $chromosome, $analysis, $cmd, $state_href ) = @_;

    print "chr=$chromosome analysis=$analysis cmd=TRUNCATE $cmd\n" if ($verbose);

    my ( $gff1_name, $delta1, $delta2, $gff2_name ) = ( $cmd =~ /(\S+)\s+([+-]*\d+)\s+([+-]*\d+)\s*=>\s*(\S+)/i );
    my $gff1_array = &get_gff( $chromosome, $gff1_name );
    my @new_gff = @{$gff1_array};

    foreach my $feature (@new_gff) {    # id, $chrom_start, $chrom_end, $chrom_strand
        if ( $feature->[3] eq '+' ) {
            $feature->[2] = $feature->[1] + $delta2;    # the end is a bit after the start
            $feature->[1] += $delta1;     # adjust the start
        }
        else {
            $feature->[1] = $feature->[2] - $delta2;    # the end is a bit after the start
            $feature->[2] -= $delta1;     # adjust the start
        }
    }

    @new_gff = sort { $a->[1] <=> $b->[1] or $a->[2] <=> $b->[2] } @new_gff;

    # store it
    $GFF_data{$gff2_name} = \@new_gff;

    print "Have ", scalar @new_gff, " FIRST_HSPs\n" if ($verbose);
}

##########################################
# delete one or more GFF lists
sub drop {
    my ( $chromosome, $analysis, $cmd, $state_href ) = @_;

    print "chr=$chromosome analysis=$analysis cmd=DROP $cmd\n" if ($verbose);
    my @gff_names = split /[,\s]+/, $cmd;
    foreach my $gff1_name (@gff_names) {
        delete $GFF_data{$gff1_name};
        print "GFF list $gff1_name has been dropped\n" if ($verbose);
    }

}

##########################################
# restrict the features to those with an ID matching a regex pattern
sub id_pattern {
    my ( $chromosome, $analysis, $cmd, $state_href ) = @_;

    print "chr=$chromosome analysis=$analysis cmd=ID_PATTERN $cmd\n" if ($verbose);

    my ( $gff1_name, $pattern, $gff2_name ) = ( $cmd =~ /(\S+)\s+(\S+)\s*=>\s*(\S+)/i );
    my $gff1_array = &get_gff( $chromosome, $gff1_name );
    my @new_gff;

    foreach my $feature ( @{$gff1_array} ) {    # id, $chrom_start, $chrom_end, $chrom_strand
        my $id = $feature->[0];
        if ( $id =~ /$pattern/ ) {
            push @new_gff, $feature;
        }
    }

    @new_gff = sort { $a->[1] <=> $b->[1] or $a->[2] <=> $b->[2] } @new_gff;

    # store it
    $GFF_data{$gff2_name} = \@new_gff;

    print "Have ", scalar @new_gff, " ID_PATTERN\n" if ($verbose);
}

##########################################
# get confirmed introns not matching a gene model
# for use when there are no pre-computed confirmed intron data

sub confirmed_introns {
    my ( $chromosome, $analysis, $cmd, $state_href ) = @_;

    print "chr=$chromosome analysis=$analysis cmd=CONFIRMED_INTRONS $cmd\n" if ($verbose);

    my ( $gff1_name, $gff2_name ) = ( $cmd =~ /(\S+)\s*=>\s*(\S+)/i );
    my $gff1_array = &get_gff( $chromosome, $gff1_name );
    my @new_gff;

    # change the exons into introns
    my @gff1_introns = $ovlp->get_intron_from_exons( @{$gff1_array} );

    # want to check the introns look reasonable
    # length > 25 bases
    # contiguous in the EST
    # ends of the intron sequence:
    # ((start eq 'gt' || start eq 'gc') && end eq 'ag') || # forward sense
    #  (start eq 'ct' && (end  eq 'ac' ||  end eq 'gc'))   # reverse sense

    my $seq = &read_chromosome($chromosome);

    foreach my $feature (@gff1_introns) {    # id, $chrom_start, $chrom_end, $chrom_strand

	my $have_homology;
	if (defined $feature->[4]) {
	    $have_homology = (abs( $feature->[5] - $feature->[4]) == 1);
	} else {
	    $have_homology = 1;
	}
        if ( ($feature->[2] - $feature->[1] > 25) &&
	     $have_homology)
        {

            # check for canonical splice sequences
            my $splice5 = uc( substr( $seq, $feature->[1] - 1, 2 ) );
            my $splice3 = uc( substr( $seq, $feature->[2] - 2, 2 ) );

#print "Feature data: $feature->[0] $feature->[1] $feature->[3] $splice5 .. $splice3\n";
#print "region: " . uc(substr($seq, $feature->[1]-3, 8)) ." .. " . uc(substr($seq, $feature->[2]-4, 8)) ."\n";
#      if ($feature->[3] eq '+') {
#	if (($splice5 ne 'GT' && $splice5 eq 'GC') || $splice3 ne 'AG') {next}
#      } else {
#	# reverse sense so $splice5 is really the 3' splice site and vice versa
#	if (($splice3 ne 'AC' && $splice3 eq 'GC') || $splice5 ne 'CT') {next}
#      }

# we are not sure of the orientation of the ESTs, so here we accept splice sites in either orientation
            if (
                (
                    ( $splice5 ne 'GT' && $splice5 eq 'GC' ) || $splice3 ne 'AG'
                )
                && ( ( $splice3 ne 'AC' && $splice3 eq 'GC' )
                    || $splice5 ne 'CT' )
              )
            {
                next;
            }

            #print "confirmed as a feature\n\n";
            push @new_gff, $feature;
        }
    }

    @new_gff = sort { $a->[1] <=> $b->[1] or $a->[2] <=> $b->[2] } @new_gff;

    # store it
    $GFF_data{$gff2_name} = \@new_gff;

    print "Have ", scalar @new_gff, " CONFIRMED_INTRONS\n" if ($verbose);

}

##########################################
# Split each feature into a 5' and a 3' end feature.
# Make new features with positions either end of the splice site at
# the start and end of the input feature.

sub split_ends {
    my ( $chromosome, $analysis, $cmd, $state_href ) = @_;

    print "chr=$chromosome analysis=$analysis cmd=SPLIT_ENDS $cmd\n" if ($verbose);

    #print "Existing GFF lists:";
    #foreach my $key (keys %GFF_data) {
    #  print "\t'$key'\n";
    #}
    #print "\n";

    my ( $gff1_name, $gff2_name, $gff3_name ) = ( $cmd =~ /(\S+)\s*=>\s*(\S+?)[,\s]+(\S+)/i );
    my $gff1_array = &get_gff( $chromosome, $gff1_name );
    my @new_5_gff;
    my @new_3_gff;

    foreach my $feature ( @{$gff1_array} ) {    # id, $chrom_start, $chrom_end, $chrom_strand
        if ( $feature->[3] eq '+' ) {
            my $start = $feature->[1];
            my $end   = $feature->[2];
            $feature->[1] = $start - 1;
            $feature->[2] = $start;
            push @new_5_gff, $feature;
            $feature->[1] = $end;
            $feature->[2] = $end + 1;
            push @new_3_gff, $feature;
        }
        else {    # '-' sense
            my $start = $feature->[1];
            my $end   = $feature->[2];
            $feature->[1] = $start - 1;
            $feature->[2] = $start;
            push @new_3_gff, $feature;
            $feature->[1] = $end;
            $feature->[2] = $end + 1;
            push @new_5_gff, $feature;
        }
    }

    @new_5_gff = sort { $a->[1] <=> $b->[1] or $a->[2] <=> $b->[2] } @new_5_gff;
    @new_3_gff = sort { $a->[1] <=> $b->[1] or $a->[2] <=> $b->[2] } @new_3_gff;

    # store them
    $GFF_data{$gff2_name} = \@new_5_gff;
    $GFF_data{$gff3_name} = \@new_3_gff;

    print "Have ", scalar @new_5_gff, " AND ", scalar @new_3_gff,
      " SPLIT_ENDS\n"
      if ($verbose);
}

##########################################
# partition splice sites depending on whether they are above or below a score threshold
sub splice_score {
    my ( $chromosome, $analysis, $cmd, $state_href ) = @_;

    print "chr=$chromosome analysis=$analysis cmd=SPLICE_SCORE $cmd\n" if ($verbose);

    my ( $gff1_name, $type, $threshold, $gff2_name, $gff3_name ) = ( $cmd =~ /(\S+)\s+(\S+)\s+([\d\.]+)\s*=>\s*(\S+?)[,\s]+(\S+)/i );
    my $gff1_array = &get_gff( $chromosome, $gff1_name );
    my @weak_gff;
    my @strong_gff;

    my $pwm = PWM->new;
    my $seq = &read_chromosome($chromosome);

    my $splice_cutoff = 0.1;

    foreach my $feature ( @{$gff1_array} ) {    # id, $chrom_start, $chrom_end, $chrom_strand
        my $strand      = $feature->[3];
        my $feature_pos = $feature->[1];

        my $score;    # score of this splice site
        if ( $strand eq '+' ) {

#$desc= "seq=". substr($seq, $feature_pos-10, 9)."/".substr($seq, $feature_pos-1, 10).  "\n";
            if ( $type eq "5'" ) {    # is it a 5' end of the intron
                $score = $pwm->splice3( $seq, $feature_pos - 1, '+' );
            }
            else {                    # or the 3' end of the intron
                $score = $pwm->splice3( $seq, $feature_pos - 1, '+' );
            }

        }
        else {                        # strand eq '-'
            if ( $type eq "5'" ) {    # is it a 5' end of the intron
                $score = $pwm->splice5( $seq, $feature_pos, '-' );
            }
            else {                    # or the 3' end of the intron
                $score = $pwm->splice5( $seq, $feature_pos, '-' );
            }
        }

        if ( $score < $splice_cutoff ) {
            push @weak_gff, $feature;
        }
        else {
            push @strong_gff, $feature;
        }

    }

    @weak_gff = sort { $a->[1] <=> $b->[1] or $a->[2] <=> $b->[2] } @weak_gff;
    @strong_gff = sort { $a->[1] <=> $b->[1] or $a->[2] <=> $b->[2] } @strong_gff;

    # store them
    $GFF_data{$gff2_name} = \@weak_gff;
    $GFF_data{$gff3_name} = \@strong_gff;

    print "Have ", scalar @weak_gff,   " WEAK SPLICE_SCORE\n"   if ($verbose);
    print "Have ", scalar @strong_gff, " STRONG SPLICE_SCORE\n" if ($verbose);
}

##########################################
# detect blastx homologies that look like frameshifts
sub frameshift {
    my ( $chromosome, $analysis, $cmd, $state_href ) = @_;

    print "chr=$chromosome analysis=$analysis cmd=FRAMESHIFT $cmd\n" if ($verbose);

    my ( $gff1_name, $gff2_name ) = ( $cmd =~ /(\S+)\s*=>\s*(\S+)/i );
    my $gff1_array = &get_gff( $chromosome, $gff1_name );
    my @frameshifted_gff;

    ###########################################
    # find the frameshifts
    ###########################################

    my $SCORE_THRESHOLD = 40;
    my $MAX_CHROM_DIFF  = 15 ; # maximum distance that the ends of the alignment are allowed to be on the chromosome
    my $MIN_CHROM_DIFF = -25 ; # minimum distance that the ends of the alignment are allowed to be on the chromosome
    my $MAX_PROT_DIFF = 15 ; # maximum distance that the ends of the HSPs are allowed to be from each other in the protein

# sort the homologies grouped by protein ID and then score and then chromosomal position
    my @homologies = sort {
             $a->[0] cmp $b->[0]
          or $a->[6] <=> $b->[6]
          or $a->[1] <=> $b->[1]
          or $a->[2] <=> $b->[2]
    } @{$gff1_array};

    # for each protein, compare all its HSPs with all its downstream HSPs
    for ( my $i = 0 ; $i < scalar @homologies ; $i++ ) {
        my $homology           = $homologies[$i];
        my $prev_protein_id    = $homology->[0];
        my $prev_chrom_start   = $homology->[1];
        my $prev_chrom_end     = $homology->[2];
        my $prev_hit_end       = $homology->[5];
        my $prev_protein_score = $homology->[6];
        if ( $prev_protein_score < $SCORE_THRESHOLD ) {
            next;
        }    # only look at high-scoring proteins

        # get the next downstream HSP
        for ( my $j = $i + 1 ; $j < scalar @homologies ; $j++ ) {
            if ( $homologies[$j]->[0] ne $prev_protein_id ) {
                last;
            }    # no longer looking at HSPs of the same protein
            if ( $homologies[$j]->[6] != $prev_protein_score ) {
                last;
            } # no longer looking at HSPs of an alignment with all the same score
            if ( $homologies[$j]->[6] < $SCORE_THRESHOLD ) {
                next;
            }    # only look at high-scoring proteins
            my $protein_id    = $homologies[$j]->[0];
            my $chrom_start   = $homologies[$j]->[1];
            my $chrom_strand  = $homologies[$j]->[3];
            my $hit_start     = $homologies[$j]->[4];
            my $protein_score = $homologies[$j]->[6];

            my $chrom_diff = $chrom_start - $prev_chrom_end;
            my $prot_diff  = $hit_start - $prev_hit_end;

            # get the difference in the frame aligned to
            # if diff % 2 is not 1 then it is a frameshift
            my $frameshift = ( abs($chrom_diff) % 3 ) - 1;

# if we are getting too far away from the end of $homologies[$i], then increment $i to stop searching it
            if ( $chrom_diff > $MAX_CHROM_DIFF ) { last; }

#print "$prev_protein_id $prev_chrom_start..$prev_chrom_end ($i) $protein_id $chrom_start.. ($j) chrom_diff $chrom_diff prot_diff $prot_diff frameshift $frameshift protein_score $protein_score prev_protein_score $prev_protein_score\n" if ($verbose);

# output any frameshifts found
# want the frame to have changed, so look at the chrom_start to prev_chrom_end difference mod 3
            if (   $frameshift
                && $chrom_diff > $MIN_CHROM_DIFF
                && $chrom_diff < $MAX_CHROM_DIFF
                && abs($prot_diff) < $MAX_PROT_DIFF )
            {

                # make a dummy GFF line for this frameshift
                my $score = ( $prev_protein_score + $protein_score ) / 2;
                my $homol;
                if ( $prev_chrom_end < $chrom_start ) {
                    $homol = [
                        $prev_protein_id, $prev_chrom_end,
                        $chrom_start,     $chrom_strand,
                        0,                0,
                        $score
                    ];
                }
                else {
                    $homol = [
                        $prev_protein_id, $chrom_start,
                        $prev_chrom_end,  $chrom_strand,
                        0,                0,
                        $score
                    ];
                }
                push @frameshifted_gff, $homol;
            }
        }
    }

    @frameshifted_gff =
      sort { $a->[1] <=> $b->[1] or $a->[2] <=> $b->[2] } @frameshifted_gff;

    # store them
    $GFF_data{$gff2_name} = \@frameshifted_gff;

    print "Have ", scalar @frameshifted_gff, " FRAMESHIFTED\n" if ($verbose);
}

##########################################
# find genes which should be merged/split based on the protein alignment

sub protein_align {
    my ( $chromosome, $analysis, $cmd, $state_href ) = @_;

    print "chr=$chromosome analysis=$analysis cmd=PROTEIN_ALIGN $cmd\n" if ($verbose);

    my ( $gff1_name, $gff2_name, $gff3_name, $gff4_name ) = ( $cmd =~ /(\S+)\s*,\s*(\S+)\s*=>\s*(\S+)\s*,\s*(\S+)/i );
    my $blastx_array = &get_gff( $chromosome, $gff1_name );
    my $CDS_array    = &get_gff( $chromosome, $gff2_name );

    my @merged_gff;    # resulting list of genes to be merged
    my @split_gff;     # resulting list of genes to be split

    #
    # find the coding exons that each protein BLASTX HSP overlaps
    #

    my @homologies = @{$blastx_array};

    my @matched = (); # the resulting list of hashes of homologies which match a coding exon
    my @matching_exons;
    my %sorted_proteins ; # protein hash keyed by protein name holding the homology and match data to find unmatched proteins that have other local matches

    my $SMALL_OVERLAP = -20;    # amount of small overlap to ignore

    my $match_obj = $ovlp->compare( $CDS_array, near => $SMALL_OVERLAP, same_sense => 0 );

    foreach my $homology (@homologies) {    # $protein_id, $chrom_start, $chrom_end, $chrom_strand

        my $protein_name  = $homology->[0];
        my $protein_score = $homology->[6];

        my $matching_exon = "";    # the name of the exon that matches;

        if ( @matching_exons = $match_obj->match($homology) ) {
            $matching_exon = $matching_exons[0][0];
            if ( grep /1/, $match_obj->matching_sense ) {

# store the homology object in a hash keyed by its name together with its 'match' status
#
# here we want a NEW COPY of the data in $homology, not just a copy of
# the reference pointing to the same data, so we force this by getting
# the array and making a new [] reference to it, then we won't have a
# problem when we push $matching_exon on the end of $homology in the
# next section
#
                my $tmp_homology = [ @{$homology} ];
                push @{$tmp_homology}, $matching_exon;
                push @{ $sorted_proteins{$protein_name} }, $tmp_homology;

                # add to the list of homologies to a coding regions
                push @{$homology}, $matching_exon;    # make a note of the exon we matched
                push @matched, $homology;
            }
        }
    }

#
# find the proteins that have patterns of similarity indicating genes should be merged or split
#

    my $prev_protein_id = ""
      ; # use to collect alignments for a protein's homology when looking for frameshifts
    my $prev_exon         = "";
    my $prev_hit_start    = -1;
    my $prev_hit_end      = -1;
    my $prev_chrom_start  = -1;
    my $prev_chrom_end    = -1;
    my $prev_score        = 0;
    my $prev_chrom_strand = "";

    my $got_a_new_exon;
    my $got_a_big_decrease_in_HSP_start;
    my $got_a_continuation_of_the_HSPs;
    my $average_blast_score;
    my $anomaly_score;

# sort the homologies grouped by protein ID and then score and then chromosomal position
    my @matched2 = sort {
             $a->[0] cmp $b->[0]
          or $a->[6] <=> $b->[6]
          or $a->[1] <=> $b->[1]
          or $a->[2] <=> $b->[2]
    } @matched;

    foreach my $homology (@matched2)
    { # $protein_id, $chrom_start, $chrom_end, $chrom_strand, $hit_start, $hit_end, $protein_score, $matching_exon, $matching_sense

        my $protein_id    = $homology->[0];
        my $chrom_start   = $homology->[1];
        my $chrom_end     = $homology->[2];
        my $chrom_strand  = $homology->[3];
        my $hit_start     = $homology->[4];
        my $hit_end       = $homology->[5];
        my $protein_score = $homology->[6];

# $homology->[7] can hold the the other_data field from the GFF line, but is usually undef
        my $matching_exon = $homology->[8];

        #  check if these are isoforms of the same gene and if so then
        #  treat them as the same gene by removing the isoform letter at
        #  the end
        if ( &isa_cds_isoform($matching_exon) ) {
            $matching_exon = &convert_cds_isoform_to_cds_id($matching_exon);
        }

       # and change any transcript names like H25P06.1.1, H25P06.1.2 to H25P06.1
        if ( &isa_transcript_isoform($matching_exon) ) {
            $matching_exon =
              &convert_transcript_isoform_to_cds_id($matching_exon);
        }

#print "Matched: $protein_id, $chrom_start, $chrom_end, $chrom_strand, $hit_start, $hit_end, $protein_score, $matching_exon\n" if ($verbose);

# look for homologies which cover two genes or genes which cover a repeated homology
        if (   $protein_id eq $prev_protein_id
            && $protein_score == $prev_score
            && $chrom_strand eq $prev_chrom_strand )
        {

            # see if exon has changed
            if ( $prev_exon ne "" && $prev_exon ne $matching_exon ) {
                $got_a_new_exon = 1;

            #print "Merge? prev $prev_exon this $matching_exon\n" if ($verbose);
            }
            else {
                $got_a_new_exon = 0;
            }

# see if the protein HSP order is unchanged - a possible case for merging if continued over into the next gene
# same sense, non-overlapping on the chromosome and jump in HSP start
#
# We want to avoid picking up the case where there are two duplicated single-exon genes
# and the homologies have a hit_start that differs by only a few bases - this should be ignored.
# So we look for cases where the start continues as greater than
# the previous mid point for evidence of a possible merge.
#
            my $prev_midpoint =
              $prev_hit_start + ( $prev_hit_end - $prev_hit_start ) / 2;

            #my $midpoint = $hit_start + ($hit_end - $hit_start)/2;
            if (
                (
                       $chrom_strand eq '+'
                    && $prev_chrom_end < $chrom_start
                    && $hit_start > $prev_midpoint
                )
                || (   $chrom_strand eq '-'
                    && $prev_chrom_end < $chrom_start
                    && $hit_start < $prev_midpoint
                    && $prev_hit_start != -1 )
              )
            {
                $got_a_continuation_of_the_HSPs = 1;
            }
            else {
                $got_a_continuation_of_the_HSPs = 0;
            }

# see if we have a big decrease in the start of the HSP - a possible case for splitting
            if (
                (
                       $chrom_strand eq '+'
                    && $prev_chrom_end < $chrom_start
                    && ( $prev_hit_start > $hit_start + 100
                        || ( $prev_hit_end > 100 && $hit_start < 30 ) )
                )
                || (
                       $chrom_strand eq '-'
                    && $prev_chrom_end < $chrom_start
                    && (
                           $prev_hit_start < $hit_start - 100
                        && $prev_hit_start != -1
                        || (   $prev_hit_start < 30
                            && $hit_start > 100
                            && $prev_hit_start != -1 )
                    )
                )
              )
            {
                $got_a_big_decrease_in_HSP_start = 1;

#print "Split? $protein_id chrom_start $chrom_end prev_chrom_start $prev_chrom_start $chrom_strand prev_hit $prev_hit_start .. $prev_hit_end hit $hit_start .. $hit_end\n" if ($verbose);
            }
            else {
                $got_a_big_decrease_in_HSP_start = 0;
            }

            # we have a merge if the pattern of HSPs shows a continuation of
            # the protein matching in order and this is a new gene and the
            # distance between the alignments is less than 5 kb
            if (   $got_a_new_exon
                && $got_a_continuation_of_the_HSPs
                && $chrom_start - $prev_chrom_end < 5000 )
            {
                push @merged_gff, $homology;
            }

            # do we have a split?
            if ( $got_a_big_decrease_in_HSP_start && !$got_a_new_exon ) {
                push @split_gff, $homology;
            }
        }

        $prev_exon         = $matching_exon;
        $prev_hit_start    = $hit_start;
        $prev_hit_end      = $hit_end;
        $prev_score        = $protein_score;
        $prev_chrom_start  = $chrom_start;
        $prev_chrom_end    = $chrom_end;
        $prev_protein_id   = $protein_id;
        $prev_chrom_strand = $chrom_strand;
    }

    @merged_gff = sort { $a->[1] <=> $b->[1] or $a->[2] <=> $b->[2] } @merged_gff;
    @split_gff = sort { $a->[1] <=> $b->[1] or $a->[2] <=> $b->[2] } @split_gff;

    # store them
    $GFF_data{$gff3_name} = \@merged_gff;
    $GFF_data{$gff4_name} = \@split_gff;

    print "Have ", scalar @merged_gff, " PROTEIN_MERGED\n" if ($verbose);
    print "Have ", scalar @split_gff,  " PROTEIN_SPLIT\n"  if ($verbose);
}
##########################################
##########################################
##########################################

# partition a list into those features which have an ID that appears
# to be from an isoform and those that don't. This is very dependent on
# the ID having enough information for you to do this and even then,
# you may have to change this routine if you are not part of the
# WormBase group where the isoform IDs have a letter appended to the
# normal sequence or CDS ID.

sub isoform {

    my ( $chromosome, $analysis, $cmd, $state_href ) = @_;

    print "chr=$chromosome analysis=$analysis cmd=ISOFORM $cmd\n" if ($verbose);

    my ( $gff1_name, $gff2_name, $gff3_name ) = ( $cmd =~ /(\S+)\s*=>\s*(\S+),\s*(\S+)/i );
    my $gff1_array = &get_gff( $chromosome, $gff1_name );
    my @isoform_list;
    my @not_isoform_list;

    foreach my $feature ( @{$gff1_array} ) {    # $id, $chrom_start, $chrom_end, $chrom_strand

        # output multiple UTR exons to the database
        my $id = $feature->[0];

        if ( &isa_cds_isoform($id) || &isa_transcript_isoform($id) ) {
            push @isoform_list, $feature;
        }
        else {
            push @not_isoform_list, $feature;
        }

    }

    # store the results
    $GFF_data{$gff2_name} = \@isoform_list;
    $GFF_data{$gff3_name} = \@not_isoform_list;

    print "Have ", scalar @isoform_list, " isoforms and ",
      scalar @not_isoform_list, " non-isoforms\n"
      if ($verbose);

}

##########################################

sub count_id {

    my ( $chromosome, $analysis, $cmd, $state_href ) = @_;

    print "chr=$chromosome analysis=$analysis cmd=COUNT_ID $cmd\n" if ($verbose);

    my ( $gff1_name, $count, $gff2_name, $gff3_name ) = ( $cmd =~ /(\S+)\s+(\d+)\s*=>\s*(\S+),\s*(\S+)/i );
    my $gff1_array = &get_gff( $chromosome, $gff1_name );
    my @lower_list;
    my @higher_list;
    my @tmp_list;
    my $prev_id = "";

    foreach my $feature ( @{$gff1_array} ) {    # $id, $chrom_start, $chrom_end, $chrom_strand

        # output multiple UTR exons to the database
        my $id = $feature->[0];

        if ( $id ne $prev_id ) {
            if ( scalar @tmp_list <= $count ) {
                push @lower_list, $feature;
            }
            else {
                push @higher_list, $feature;
            }
            @tmp_list = ();
        }

    }

    # store the results
    $GFF_data{$gff2_name} = \@lower_list;
    $GFF_data{$gff3_name} = \@higher_list;

    print "Have ", scalar @lower_list, " <= $count and ", scalar @higher_list,
      " > $count\n"
      if ($verbose);

}

##########################################

sub read_gff {

    my ( $chromosome, $analysis, $cmd, $state_href ) = @_;

    print "chr=$chromosome analysis=$analysis cmd=READ_GFF $cmd\n" if ($verbose);

    my ( $file_name, $source, $type, $id_after, $gff_name ) = ( $cmd =~ /(\S+)\s+(\S+)\s+(\S+)\s+(\S+)\s*=>\s*(\S+)/i );

    $file_name =~ s/\<CHROMOSOME\>/$chromosome/g;

#print "Reading GFF file: $file_name source: $source, type: $type, id_after: $id_after\n";

    my %GFF_file_data = (
        path       => $file_name,    # use this explicit path for the file
        gff_source => $source,
        gff_type   => $type,
        ID_after   => $id_after,
    );

    my @gff = $ovlp->read_GFF_file( \%GFF_file_data, $chromosome );

    # store the results
    $GFF_data{$gff_name} = \@gff;

    print "Have read " . scalar @gff . " lines for '$gff_name'\n" if ($verbose);

}

##########################################
# get a GFF list if it has not already been got
# stores the GFF list in %GFF_data and also returns an aref to the list
sub get_gff {
    my ( $chromosome, $gff_name ) = @_;

    print "Getting GFF list '$gff_name' for chromosome $chromosome\n" if ($verbose);

    if ( !exists $GFF_data{$gff_name} ) {
        print "Reading GFF data for '$gff_name'\n" if ($verbose);
        my @gff = $ovlp->get( $gff_name, $chromosome );
        $GFF_data{$gff_name} = \@gff;
        print "Have read " . scalar @gff . " lines for '$gff_name'\n" if ($verbose);
    }

    return $GFF_data{$gff_name};
}
##########################################
##########################################
##########################################
##########################################
##########################################
##########################################
##########################################
##########################################
##########################################
##########################################
# read the chromosmal sequence
#  my $seq = read_chromosome($chromosome);

sub read_chromosome {
    my $chromosome = shift;

    # if we have already read in the sequence entries, return the one for this chromosome
    if ( exists $dna_entry{$chromosome} ) { return $dna_entry{$chromosome}; }

    my ($seq_file) = glob("$dnadir/$chromosome.dna");

    my $seq = &read_file($seq_file);

    if ( !defined $seq ) {
        $seq = &read_entry( glob("$dnadir"), $chromosome );
    }

    return $seq;

}

##########################################
# read file

sub read_file {
    my ($file) = @_;

    my $seq;

    # try to open a single-chromosome entry file e.g. for elegans, briggsae
    if ( open( SEQ, $file ) ) {
        my $tmp = $/;
        $/   = "";
        $seq = <SEQ>;
        close SEQ;
        $/ = $tmp;

        $seq =~ s/>.*\n//;    # remove one title line
        $seq =~ s/\n//g;
    }
    else {
        print "Can't open the dna file for $file : $!\n" if ($verbose);
    }

    return $seq;
}

##########################################
# read entry in a file if the file consists of lots of entries

sub read_entry {
    my ( $dir, $chromosome ) = @_;

# if we have already read in the sequence entries, return the one for this chromosome
    if ( exists $dna_entry{$chromosome} ) { return $dna_entry{$chromosome}; }

    my $file = "$dir/supercontigs.fa";    # it might be called this
    print "Trying: $file\n" if ($verbose);
    if ( !-e $file ) {
        print "Can't open the dna file for $file : $!\n" if ($verbose);
        $file = "$dir/${species}.fa"
          ;   # or it might be called this - the nomenclature is not yet certain
        print "Trying: $file\n" if ($verbose);

    }
    elsif ( !-e $file ) {
        die "Can't find chromosome sequence file for $dir, $chromosome\n";
    }

    print "Reading DNA sequence entries\n" if ($verbose);
    my $seq;

    if ( open( SEQ, $file ) ) {
        my $entry;
        while ( my $line = <SEQ> ) {
            chomp $line;
            if ( $line =~ /^>(\S+)/ ) {
                if ( defined $entry ) {    # store the previous sequence entry
                    $dna_entry{$entry} = $seq;
                }
                $seq   = "";
                $entry = $1;
            }
            else {
                $seq .= $line;
            }
        }
        if ( defined $entry ) {            # store the last entry
            $dna_entry{$entry} = $seq;
        }
        close SEQ;
    }
    else {
        die "Can't open the dna file for $file : $!\n";
    }

    if ( !exists $dna_entry{$chromosome} ) {
        die "No DNA sequence was found for $chromosome\n";
    }
    return $dna_entry{$chromosome};
}

####################################################################################
#
# General Subroutines
#
####################################################################################

####################################################################################
#
# Returns True if the input CDS ID is an isoform

sub isa_cds_isoform {

    my $CDS = shift;
    return ( $CDS =~ /(\S+\.\d+)[a-z]/ );

}

####################################################################################
#
# Converts a CDS isoform ID into the corresponding base ID of any CDS for this gene

sub convert_cds_isoform_to_cds_id {

    my $CDS = shift;
    $CDS =~ /(\S+\.\d+)[a-z]/;
    return $1;

}

####################################################################################
#
# Returns True if the input transcript ID is an isoform

sub isa_transcript_isoform {

    my $TRANS = shift;
    return ( $TRANS =~ /(\S+\.\d+)\.\d+/ );

}

####################################################################################
#
# Converts a transcript isoform ID into the corresponding base ID of any CDS for this gene

sub convert_transcript_isoform_to_cds_id {

    my $TRANS = shift;
    $TRANS =~ /(\S+\.\d+)\.\d+/;
    return $1;

}
####################################################################################
# checks if there is sense = + or - and if not, outputs two anomaly records, one for each sense.
# else it just outputs one record for the specified sense.

sub output_to_database {
    my (
        $anomaly_type, $chromosome,   $anomaly_id,    $chrom_start,
        $chrom_end,    $chrom_strand, $anomaly_score, $explanation
    ) = @_;

  # get the lab for this location
  # we don't use this is the public version, always set lab and clone to a blank
    my $lab         = '';
    my $clone       = '';
    my $clone_start = 0;
    my $clone_end   = 0;

    # calculate the window value as blocks of 10 kb
    my $window = int( $chrom_start / 10000 );

# truncate the $anomaly_id string to 32 characters (the width of the 'thing_id' field in the 'anomaly' table)
    $anomaly_id = substr( $anomaly_id, 0, 32 );

    # if there is no strand information, put one anomaly in for each strand
    if ( $chrom_strand ne '+' && $chrom_strand ne '-' ) {
        &put_anomaly_record_in_database(
            $anomaly_type, $chromosome, $anomaly_id,  $chrom_start,
            $chrom_end,    '+',         $window,      $anomaly_score,
            $explanation,  $clone,      $clone_start, $clone_end,
            $lab
        );
        &put_anomaly_record_in_database(
            $anomaly_type, $chromosome, $anomaly_id,  $chrom_start,
            $chrom_end,    '-',         $window,      $anomaly_score,
            $explanation,  $clone,      $clone_start, $clone_end,
            $lab
        );

    }
    else {
        &put_anomaly_record_in_database(
            $anomaly_type, $chromosome,   $anomaly_id,  $chrom_start,
            $chrom_end,    $chrom_strand, $window,      $anomaly_score,
            $explanation,  $clone,        $clone_start, $clone_end,
            $lab
        );

    }
}

##########################################
# output the record to the data file.

sub put_anomaly_record_in_database {

    my (
        $anomaly_type, $chromosome,   $anomaly_id,  $chrom_start,
        $chrom_end,    $chrom_strand, $window,      $anomaly_score,
        $explanation,  $clone,        $clone_start, $clone_end,
        $lab
    ) = @_;

    # ignore this if the score is less than 0.01
    if ( $anomaly_score < 0.01 ) { return; }

    #################################
    # output the data to the GFF file
    #################################

#  if ($supplementary) {
#    print OUTPUT_GFF "$chromosome\tcuration_anomaly\t$anomaly_type\t$chrom_start\t$chrom_end\t$anomaly_score\t$chrom_strand\t.\tEvidence \"$anomaly_id\"\n";
#  }

    #############################################
    # write the anomaly data to the data file
    #############################################

    # count the anomalies found
    $anomaly_count{$anomaly_type}++;

    # and output the line for the datafile
    print DAT "INSERT\t$anomaly_type\t$chromosome\t$anomaly_id\t$chrom_start\t$chrom_end\t$chrom_strand\t$anomaly_score\t$explanation\t$clone\t$clone_start\t$clone_end\t$lab\n";

}

##########################################
# Delete things that have not been updated in this run that you
# would expect to have been updated like protein-homology-based
# anomalies that might have gone away.  This also deletes anomalies
# that we are no longer putting into the database and which can now be
# removed.
#
# This means that an anomaly based on something like a similarity to a
# protein where the protein is no longer existing will be removed from
# the database

sub delete_anomalies {

    my ($type) = @_;

    # output the line for the datafile
    print DAT "DELETE\t$type\n";

}

##########################################

sub usage {
    my $error = shift;

    if ( $error eq "Help" ) {

        # Normal help menu
        system( 'perldoc', $0 );
        exit(0);
    }
}

##########################################

# Add perl documentation in POD format
# This should expand on your brief description above and
# add details of any options that can be used with the program.
# Such documentation can be viewed using the perldoc command.

__END__


=pod

=head2 NAME - find_anomalies.pl

=head1 USAGE

=over 4

=item find_anomalies.pl  [-options]

=back

This writes a data file ready to populate the worm_anomaly mysql
database with data describing some types of anomalies.

It can be run periodicaly e.g. every build which is especially useful
if there have been corrections made to the genomic sequence

These anomalies can be inspected using the script:

history_maker.pl -anomalies -chromosome X

script_template.pl MANDATORY arguments:

=over 4

=item none

=back

script_template.pl  OPTIONAL arguments:

=over 4

=item -species

The species - the default is 'elegans'

=over 4

=item -h, Help

=back

=over 4
 
=item -debug, Verbose/Debug mode
 
=back

=over 4

=item -test, Test mode

=back

=over 4
    
=item -verbose, output lots of chatty test messages

=back
                                                                                             

=head1 REQUIREMENTS

=over 4

=item GFF files for the input data

=back

=head1 AUTHOR

=over 4

=item Gary Williams (gw3@sanger.ac.uk)

=back

=cut


