#
# Perl Module for GFF
#
# Cared for by Tim Hubbard <th@sanger.ac.uk>
#
# Copyright Genome Research Limited (1997). Please see information on licensing in LICENSE

package GFF;

use vars qw($AUTOLOAD @ISA @EXPORT_OK);
use Exporter;
use Carp;
use strict;
use GeneFeature;
use HomolGeneFeature;
use GeneFeaturePair;
use GFFPair;

#
# Place functions/variables you want to *export*, ie be visible from the caller package into @EXPORT_OK
#

@EXPORT_OK = qw();

#
# @ISA has our inheritance.
#

@ISA = ( 'Exporter' );

my %fields = (
	      feature => undef, # going to be an array of GeneFeatures
	      filename => undef
	      );

#
# creates a new object
#

sub new {
    my $ref = shift;
    my $class = ref($ref) || $ref;
    
    my $self = {
	'_permitted' => \%fields,
	%fields, };

    $self->{'feature'} = []; # makes feature an anonymous array.
    bless $self, $class;
    return $self;
}

#
# Autoload up my field name set/get
#

sub AUTOLOAD {
    my $self = shift;
    my $type = ref($self) || carp "Ok, $self is not an object of mine!";
    my $name = $AUTOLOAD;

    # don't propagate DESTROY messages...

    $name =~ /::DESTROY/ && return;

    $name =~ s/.*://; #get only the bit we want
    unless (exists $self->{'_permitted'}->{$name} ) {
	carp "In type $type, can't access $name - probably passed a wrong variable into Sequence  ";
    }
    if (@_) {
	return $self->{$name} = shift;
    } else {
	return $self->{$name}
    }
}

#
# Method to add a reference to a GeneFeature object to a GFF object
#

sub addGeneFeature {
    my $self = shift;
    my $gf = shift;
    my $function = shift;

    if(!$function || &$function($gf)){
	push( @{$self->{'feature'}}, $gf);
    }
}

#
# Method to append GeneFeatures in a GFF object to another GFF object
#

sub addGFF {
    my $self = shift;
    my $other = shift;
    my $function = shift;

    my($gf);

    foreach $gf ($other->eachGeneFeature()){
	if(!$function || &$function($gf)){
	    push( @{$self->{'feature'}}, $gf);
	}
    }
}

# 
# Method to return an array of refs to the GeneFeatures in a GFF object
#

sub eachGeneFeature {
    my $self = shift;
    my @array;
    my $gf;

    foreach $gf ( @{$self->{'feature'}} ){
	push(@array,$gf);
    }
    return @array;
}

#
# Method to create GeneFeature object for each line of a stream from a GFF format file
# and add references to them to the GFF object
#

sub read {
    my $self = shift;
    my $file = shift;
    my $function = shift;
    my $function2 = shift;
    my ($gf,$count);

    $count = 0;

    if( !defined $file ) {
	carp("sorry - can't read a GFF with no file!");
	return 0;
    }

    # loop over file, make new gf objects and add them

    while( <$file> ) {
	chomp($_);

	/^#/ && next;

# convert file format using function
	if($function2){
	    $_=&$function2($_);
	}

	$gf = &GeneFeature::new_from_line($_);

	$self->addGeneFeature($gf,$function);
	
	$count++;
    }

    return $count;
}

#
# Method to create GeneFeature object for each line of a stream from an MSP file
# (output from MSPcrunch) and add references to them to the GFF object
#

sub read_msp {
    my $self = shift;
    my $file = shift;
    my $source = shift;
    my $function = shift;
    my $function2 = shift;
    my ($gf,$count);

    if( !defined $file ) {
	carp("sorry - can't read a GFF with no file!");
	return 0;
    }

    # loop over file, make new gf objects and add them
    $count = 0;
    while( <$file> ) {
	chomp($_);
	/^#/ && next;
	$gf = &HomolGeneFeature::new_from_msp($_,$source,$function2);
	$self->addGeneFeature($gf,$function);
	$count++;
    }
    return $count;
}

#
# Method to dump out a GFF object, via method to dump out a GeneFeature object.
#
    
sub dump {
    my $self = shift;
    my $file = shift;
    my $tab = shift;
    my $gf;

    foreach $gf ( $self->eachGeneFeature() ) {
	$gf->dump($file,$tab);
    }
    
}

#
# Method to dump out a GFF object, via method to dump out a GeneFeature object.
# Includes information about any pairing.
#
    
sub dump_pair {
    my $self = shift;
    my $file = shift;
    my $tab = shift;
    my $gf;

    foreach $gf ( $self->eachGeneFeature() ) {
	$gf->dump_pair($file,$tab);
    }
    
}

#
# Method to generate a new GFF object based on a filtered version of the object.
# Filtering is carried out by passing a reference to a subroutine, which is applied
# to each GeneFeature object.
#

sub filter {
    my $self = shift;
    my $function = shift;
    my $gf;

    my $new = new GFF;
    
    foreach $gf ($self->eachGeneFeature()){
	$new->addGeneFeature($gf,$function);
    }
    $new;
}

#
# Method to count the number of GeneFeature objects in a GFF object
#

sub count {
    my $self = shift;
    my $gf;
    my $count;

    $count=0;
    foreach $gf ($self->eachGeneFeature()){
	$count++;
    }
    $count;
}

#
# Method return strand of gf objects in gff
#

sub strand {
    my $self = shift;
    my ($strand,%strand,$strandall,$gf);

    foreach $gf ($self->eachGeneFeature()){
	$strand{$gf->strand}++;
    }
    foreach $strand (sort keys %strand){
	$strandall.=$strand;
    }
    return $strandall;
}

#
# Method to generate a new GFF object based on a filtered version of the object
# based on any pairwise overlap.  Option to report overlaps to filehandle.
#

sub self_overlap {
    my $self = shift;
    my $verbose = shift;

    my (@gf,@del);
    my ($i,$j,$gf,$string);

# setup array of GeneFeature objects
    @gf=$self->eachGeneFeature();

# loop over each
    for($i=0;$i<scalar(@gf);$i++){
	$gf=$gf[$i];

# loop over each and build list of deletes
	undef @del;
	for($j=$i+1;$j<scalar(@gf);$j++){
	    if($string=$gf->overlap_logical($gf[$j],$verbose)){
		print $verbose "  OVERLAP: $string (DELETE)\n" if $verbose;
		push(@del,$j);
	    }
	}

# apply deletes from this pass
	foreach $j (@del){
	    splice(@gf,$j,1);
	}
    }

# build new object from remaining list
    my $new = new GFF;
    foreach $gf (@gf){
	$new->addGeneFeature($gf);
    }
    $new;

}

#
# Method to generate a new GFF object based on intersection with second GFF object
# tolerance specifies error at each edge: tolerance=-1 specifies any overlap;
# tolerance=0 specifies exact match.
#

sub intersect_overlap {
    my $self = shift;
    my $other = shift;
    my $tolerance = shift;
    my $single = shift;
    my $strand = shift;

    my (@gf,@string);
    my ($gf,$gf2,$string,%saved);

# new object from intersection
    my $new = new GFF;

# setup array of GeneFeature objects
# from GFF or GeneFeature
    if(ref($other) eq "GFF"){
	@gf=$other->eachGeneFeature();
    }elsif(ref($other) eq "GeneFeature"){
	@gf=($other);
    }elsif(ref($other) eq "HomolGeneFeature"){
	@gf=($other);
    }else{
	die "$other if of unexpected class";
    }

# loop over each
    foreach $gf ($self->eachGeneFeature()){
# loop over other 
	foreach $gf2 (@gf){
	    @string=$gf->match($gf2,$tolerance,$single,$strand);
	    if($string[0]){
		unless($saved{$gf}){
		    $new->addGeneFeature($gf);
		    $saved{$gf}=1;
		}
	    }
	}
    }
    return $new;
}

#
# Method to generate a new GFFPair object based on intersection with second GFF object
# tolerance specifies error at each edge: tolerance=-1 specifies any overlap;
# tolerance=0 specifies exact match.
# single=1 specifies that only one edge needs to match tolerance
# strand=1 specifies that orientation must be identical
#

sub intersect_overlap_pair {
    my $self = shift;
    my $other = shift;
    my $tolerance = shift;
    my $single = shift;
    my $strand =shift;
    my $verbose = shift;
    my $file = shift;

    my (@gf,@string);
    my ($gf,$gf1,$gf2,$string);

# new object from intersection
    my $new = &GFFPair::new('GFFPair');

# setup array of GeneFeature objects
    @gf=$other->eachGeneFeature();

# loop over each
    foreach $gf1 ($self->eachGeneFeature()){

# loop over other 
	foreach $gf2 (@gf){
	    @string=$gf1->match($gf2,$tolerance,$single,$strand);
#	    print "$string[0],$string[3]\n";
	    if($string[0]){
# create new GeneFeaturePair object
		$gf = &GeneFeaturePair::new_from_commandline($gf1,$gf2,$string[1],$string[2]);
# add to GFFPair object
		$new->addGeneFeature($gf);
# add ref to GeneFeaturePair
		$gf1->addPair($gf);
		$gf2->addPair($gf);
# report strand orientation differences
	    }elsif($strand && $string[3]==3 && $verbose){
		print $file "  ERROR: DIRECTION? ".$string[4]."\n";
	    }
	}
    }
    return $new;

}

#
# Method to generate a new GFF object which is a union of 2 GFF objects
#

sub union {
    my $self = shift;
    my $other = shift;

    my $new = new GFF;
    my $gf;

# first add everything in 1st object
    foreach $gf ($self->eachGeneFeature()){
	$new->addGeneFeature($gf);
    }

# second add everything in 2nd object but not in 1st object
    foreach $gf ($other->eachGeneFeature()){
	unless($self->member($gf)){
	    $new->addGeneFeature($gf);
	}
    }
    return $new;
}

#
# Method to generate a new GFF object which is everything in first GFF that is not in second GFF
#

sub intersect_not {
    my $self = shift;
    my $other = shift;

    my $new = new GFF;
    my $gf;

# add everything in 1st object unless it is in the 2nd object
    foreach $gf ($self->eachGeneFeature()){
	unless($other->member($gf)){
	    $new->addGeneFeature($gf);
	}
    }
    return $new;
}

#
# Method to test if a GeneFeature object is a member of an existing GFF object
# (note - comparison is by reference only - does not check if 2 GeneFeature objects
#  contain the same data!)
#

sub member {
    my $self = shift;
    my $element = shift;
    my ($gf,$flag);
    $flag=0;
    foreach $gf ($self->eachGeneFeature()){
	if($gf == $element){
	   $flag=1;
	   last;
       }
    }
    return $flag;
}

#
# Method to find max and min range of members of a GFF object.
#

sub max_min_range {
    my $self = shift;
    my ($gf,$minrng,$maxrng);
    $maxrng=0;
    $minrng=1000000000;
    foreach $gf ($self->eachGeneFeature()){
	$maxrng=&GeneFeature::max($maxrng,$gf->end);
	$minrng=&GeneFeature::min($minrng,$gf->start);
    }
    ($maxrng,$minrng);
}

#
# Method to find min range of members of a GFF object.
#

sub start_range {
    my $self = shift;
    my ($gf,$minrng);
    $minrng=1000000000;
    foreach $gf ($self->eachGeneFeature()){
	$minrng=&GeneFeature::min($minrng,$gf->start);
    }
    return $minrng;
}

#
# Method to find max range of members of a GFF object.
#

sub end_range {
    my $self = shift;
    my ($gf,$maxrng);
    $maxrng=0;
    foreach $gf ($self->eachGeneFeature()){
	$maxrng=&GeneFeature::max($maxrng,$gf->end);
    }
    return $maxrng;
}

#
# Method to find max score of members of a GFF object.  If $paired, then only evaluates
# for objects which are lined to GFFPair object
#

sub maxScore {
    my $self = shift;
    my $paired = shift;
    my ($gf,$mxscore);
    $mxscore=0;
    foreach $gf ($self->eachGeneFeature()){
	if(!$paired || $gf->getPair_logical){
	    $mxscore=&GeneFeature::max($mxscore,$gf->score);
	}
    }
    return $mxscore;
}

#
# Method to return GF object with max score in a GFF object.
#

sub maxScoreGf {
    my $self = shift;
    my $function = shift;
    my ($gf,$mxscore,$mxgf);
    $mxscore=0;
    foreach $gf ($self->eachGeneFeature()){
	if(!$function || &$function($gf)){
	    if($gf->score>$mxscore){
		$mxscore=$gf->score;
		$mxgf=$gf;
	    }
	}
    }
    return $mxgf;
}

#
# Method to build an array of GFF objects each containing a group of
# GeneFeature objects sharing some attribute [pairwise comparison]
#

sub cluster {
    my $self = shift;
    my $function = shift;
    my $single = shift;

    my (@gf);
    my ($i,$j,$gf,$gf1,$gf2,$ngf,
	$nc1,$nc2,$nc,%cluster,$ncluster,
	$cluster,@cluster);

# setup array of GeneFeature objects
    @gf=$self->eachGeneFeature();
    $ngf=scalar(@gf);
    $ncluster=0;

# loop over each
    for($i=0;$i<$ngf;$i++){
	$gf1=$gf[$i];

# loop over each other
	for($j=$i+1;$j<$ngf;$j++){
	    $gf2=$gf[$j];

# test relationship between $gf[$i] and $gf[j]
	    if(&$function($gf1,$gf2)){

# if both in existing cluster, merge clusters
		$nc1=$cluster{$gf1};
		$nc2=$cluster{$gf2};
		if($nc1 && $nc2){
# merge clusters if identical
		    if($nc1!=$nc2){
			push(@$nc1,@$nc2);
			foreach $cluster (keys %cluster){
			    if($cluster{$cluster}==$nc2){
				$cluster{$cluster}=$nc1;
			    }
			}
			undef @$nc2;
			$ncluster--;
		    }
# nc1 exists -> add nc2
		}elsif($nc1){
		    push(@$nc1,$gf2);
		    $cluster{$gf2}=$nc1;
		}elsif($nc2){
		    push(@$nc2,$gf1);
		    $cluster{$gf1}=$nc2;
		}else{
		    $nc=[];
		    push(@cluster,$nc);
		    push(@$nc,$gf1);
		    push(@$nc,$gf2);
		    $cluster{$gf1}=$nc;
		    $cluster{$gf2}=$nc;
		    $ncluster++;
		}
	    }
	}
    }

# make clusters out of each singleton
    if($single){
	foreach $gf (@gf){
	    unless($cluster{$gf}){
		$nc=[];
		push(@cluster,$nc);
		push(@$nc,$gf);
	    }
	}
    }

# build an array of new objects from set of clusters
    undef @gf;
    foreach $cluster (@cluster){
	my $new = new GFF;
	foreach $gf (@$cluster){
	    $new->addGeneFeature($gf);
	}
	push(@gf,$new);
    }
    @gf;
}

#
# Method to build an array of GFF objects each containing a group of
# GeneFeature objects sharing some attribute [fixed, named]
#

sub group {
    my $self = shift;
    my $function = shift;

    my (%gf);
    my ($gf,%group,$group,$name);

    foreach $gf ($self->eachGeneFeature()){
	if($name=&$function($gf)){
	    push(@{$group{$name}},$gf);
	}
    }

# build an array of new objects from set of clusters
    foreach $group (keys %group){
	my $new = new GFF;
	foreach $gf (@{$group{$group}}){
	    $new->addGeneFeature($gf);
	}
	$gf{$group}=$new;
    }
    return %gf;
}

#
# Method to add a label to each GeneFeature in this GFF point back to this GFF
#

sub label {
    my $self = shift;
    my $membergroup = shift;
    my ($gf);
    foreach $gf ($self->eachGeneFeature()){
	$gf->addMember($self,$membergroup);
    }
}

#
# Method to make a list of the corresponding list of GFF's that the GeneFeatures in
# this GFF are paired with
#

sub label_pair {
    my $self = shift;
    my $membergroup = shift;
    my ($gf1,$gf2,$gfp,%gff,$gff,@gff);
# loop over each gf
    foreach $gf1 ($self->eachGeneFeature()){
# loop over each gfp
	foreach $gfp ($gf1->getPairs){
# if find a pairing, follow it...
	    if($gfp){
# get paired GeneFeature
		$gf2=$gfp->other($gf1);
# look for reference back to to GFF
		$gff=$gf2->getMember($membergroup);
		if($gff){
		    unless($gff{$gff}){
			push(@gff,$gff);
		    }
		    $gff{$gff}++;
		}
	    }
	}
    }
#    return %gff;
    return @gff;
}

#
# Method to return a number of paired GeneFeature objects between 2 GFF's
#

sub shared_paired {
    my $self = shift;
    my $other = shift;
    my ($gf1,$gf2,$gfp,%gf,$gf,$n);
# loop over each gf
    foreach $gf ($other->eachGeneFeature()){
	$gf{$gf}=1;
    }
    $n=0;
    foreach $gf1 ($self->eachGeneFeature()){
# loop over each gfp
	foreach $gfp ($gf1->getPairs){
# if find a pairing, follow it...
	    if($gfp){
# get paired GeneFeature
		$gf2=$gfp->other($gf1);
		if($gf{$gf2}){
		    $n++;
		}
	    }
	}
    }
    $n;
}

#
# Method to return a hash of the cases and their occurance for a feature
#

sub features {
    my $self = shift;
    my $function = shift;
    my (%hash,$gf,$tmp);
# loop over each gf
    foreach $gf ($self->eachGeneFeature()){
	$tmp=&$function($gf);
	$hash{$tmp}++;
    }
    %hash;
}

#
# Method to test if features returns a single value and return it
#

sub theFeature {
    my $self = shift;
    my $function = shift;
    my ($feature,%hash,$n,$ref);
    $n=0;
    %hash=$self->features($function);
    foreach $feature (keys %hash){
	$ref=$feature;
	$n++;
    }
    if($n!=1){
	die "Error in GFF::theFeature - more than one feature found ($n)";
    }
    $ref,$hash{$ref};
}

1;  # says use was ok
__END__

=head1 NAME

GFF.pm

=head1 DESCRIPTION

B<GFF> is a Perl Object for Gene Finding Feature format (see
http://www.sanger.ac.uk/Users/rd/gff.html for GFF specification).  A
GFF object is a set of GeneFeature objects (see GeneFeature.pm).

=head1 AUTHOR

B<Tim Hubbard> Email th@sanger.ac.uk

=head1 Overview

=for html
<a href="GeneFeature.html">GeneFeature.pm</a> Object for a single gene feature.<br>
<a href="HomolGeneFeature.html">HomolGeneFeature.pm</a> Homology gene feature 
(object class derived from GeneFeature.pm).<br>
<a href="GFF.html">GFF.pm</a> Container object for GeneFeature or HomolGeneFeature objects.<br>
<a href="GeneFeaturePair.html">GeneFeaturePair.pm</a> Object to sort links between 
GeneFeature objects (uncertain usefulness and may be deprecated in the future)<br>
<a href="GFFPair.html">GFFPair.pm</a> Container object for GeneFeaturePair objects
(uncertain usefulness and may be deprecated in the future)<br>

=head1 Methods to read data into and output data from GFF objects.

=over

=item read

Create GeneFeature object for each line of a stream from a GFF format
file and add references to them to the GFF object

=item read_msp

Create HomolGeneFeature object for each line of a stream from an MSP
file (output from MSPcrunch) and add references to them to the GFF
object

=item dump

Dump out a GFF object (via method to dump out a GeneFeature object).

=item dump_pair

Dump out a GFF object (via method to dump out a GeneFeature object).
Includes information about any pairing.

=back

=head1 Methods to combine/manipulate GFF objects

=over

=item addGeneFeature

Add a reference to a GeneFeature object to a GFF object

=back

=head1 Methods to do calculations on GFF objects

=over

=back

=head1 Misc methods.

=over

=item eachGeneFeature

Method to return an array of refs to the GeneFeatures in a GFF object

=item filter

Method to generate a new GFF object based on a filtered version of the
object.  Filtering is carried out by passing a reference to a
subroutine, which is applied to each GeneFeature object.

=item count

Method to count the number of GeneFeature objects in a GFF object

=item self_overlap

Method to generate a new GFF object based on a filtered version of the
object based on any pairwise overlap.  Option to report overlaps to
filehandle.

=item intersect_overlap

Method to generate a new GFF object based on intersection with second
GFF object tolerance specifies error at each edge: tolerance=-1
specifies any overlap; tolerance=0 specifies exact match.

=item union

Method to generate a new GFF object which is a union of 2 GFF objects

=item member

Method to test if a GeneFeature object is a member of an existing GFF
object (note - comparison is by reference only - does not check if 2
GeneFeature objects contain the same data!)

=item addGFF

Method to append GeneFeatures in a GFF object to another GFF object

=item cluster

Method to build an array of GFF objects each containing a group of
GeneFeature objects sharing some attribute [pairwise comparison]

=item end_range

Method to find max range of members of a GFF object.

=item features

Method to return a hash of the cases and their occurance for a feature

=item group

Method to build an array of GFF objects each containing a group of
GeneFeature objects sharing some attribute [fixed, named]

=item intersect_not

Method to generate a new GFF object which is everything in first GFF that is not in second GFF

=item intersect_overlap_pair

Method to generate a new GFFPair object based on intersection with second GFF object
tolerance specifies error at each edge: tolerance=-1 specifies any overlap;
tolerance=0 specifies exact match.
single=1 specifies that only one edge needs to match tolerance
strand=1 specifies that orientation must be identical

=item label

Method to add a label to each GeneFeature in this GFF point back to this GFF

=item label_pair

Method to make a list of the corresponding list of GFF's that the GeneFeatures in
this GFF are paired with

=item maxScore

Method to find max score of members of a GFF object.  If $paired, then only evaluates
for objects which are lined to GFFPair object

=item maxScoreGf

Method to return GF object with max score in a GFF object.

=item max_min_range

Method to find max and min range of members of a GFF object.

=item shared_paired

Method to return a number of paired GeneFeature objects between 2 GFF's

=item start_range

Method to find min range of members of a GFF object.

=item strand

Method return strand of gf objects in gff

=item theFeature

Method to test if features returns a single value and return it

