#!/usr/bin/perl -w
#-----------------------------------------------------------+
#                                                           |
# Name.pl                                                   |
#                                                           |
#-----------------------------------------------------------+
#                                                           |
#  AUTHOR: James C. Estill                                  |
# CONTACT: JamesEstill_@_gmail.com                          |
# STARTED: 00/00/2007                                       |
# UPDATED: 00/00/2007                                       |
#                                                           |
# DESCRIPTION:                                              |
#  Short Program Description                                |
#                                                           |
# USAGE:                                                    |
#  ShortFasta Infile.fasta Outfile.fasta                    |
#                                                           |
# VERSION: $Rev$                                            |
#                                                           |
#                                                           |
# LICENSE:                                                  |
#  GNU General Public License, Version 3                    |
#  http://www.gnu.org/licenses/gpl.html                     |  
#                                                           |
#-----------------------------------------------------------+
# TO DO:
# Sub selection of each scaffold in the GFF file.
# This is using GFF file and object
# "Bio::SeqFeature::Annotated
# Methods documented at
# http://doc.bioperl.org/bioperl-live/Bio/SeqFeature/Annotated.html
# It would be better to load these to a GFF data store and use features
# to fetch the overlapping features. This will be more scalable
# http://search.cpan.org/~cjfields/BioPerl-1.6.1/Bio/DB/SeqFeature/Segment.pm
# 
package DAWGPAWS;

#-----------------------------+
# INCLUDES                    |
#-----------------------------+
use Graph;
use strict;
use warnings;
use Getopt::Long;
# The following needed for printing help
use Pod::Select;               # Print subsections of POD documentation
use Pod::Text;                 # Print POD doc as formatted text file
use IO::Scalar;                # For print_help subfunction
use IO::Pipe;                  # Pipe for STDIN, STDOUT for POD docs
use File::Spec;                # Convert a relative path to an abosolute path
#use Set::IntervalTree;
#use Data::Range::Compare;      #
#use Bio::Range;
use Bio::FeatureIO;

#-----------------------------+
# Program VARIABLES           |
#-----------------------------+
my ($VERSION) = q$Rev$ =~ /(\d+)/;
# Get GFF version from environment, GFF2 is DEFAULT
my $gff_ver = uc($ENV{DP_GFF}) || "GFF2";
my %all_sources = ();          # Sources

# The parent feature to fetch overlaps for, others 
# that may make use of this include features such as gene or transposable_element
# that may have child features you want to fetch. This may require
# moving to bioperl GFF feature objects.
my $parent_feature = "LTR_retrotransposon";

#-----------------------------+
# VARIABLE SCOPE              |
#-----------------------------+
my $infile;
my $outfile;

# BOOLEANS
my $quiet = 0;
my $verbose = 0;
my $show_help = 0;
my $show_usage = 0;
my $show_man = 0;
my $show_version = 0;
my $do_test = 0;                  # Run the program in test mode

#-----------------------------+
# COMMAND LINE OPTIONS        |
#-----------------------------+
my $ok = GetOptions(# REQUIRED OPTIONS
		    "i|infile=s"  => \$infile,
                    "o|outfile=s" => \$outfile,
		    "f|feature=s" => \$parent_feature,
		    # ADDITIONAL OPTIONS
		    "gff-ver=s"   => \$gff_ver,
		    "q|quiet"     => \$quiet,
		    "verbose"     => \$verbose,
		    # ADDITIONAL INFORMATION
		    "usage"       => \$show_usage,
		    "test"        => \$do_test,
		    "version"     => \$show_version,
		    "man"         => \$show_man,
		    "h|help"      => \$show_help,);

 
#-----------------------------+
# STANDARDIZE GFF VERSION     |
#-----------------------------+
unless ($gff_ver =~ "GFF3" || 
	$gff_ver =~ "GFF2") {
    # Attempt to standardize GFF format names
    if ($gff_ver =~ "3") {
	$gff_ver = "GFF3";
    }
    elsif ($gff_ver =~ "2") {
	$gff_ver = "GFF2";
    }
    else {
	print "\a";
	die "The gff-version \'$gff_ver\' is not recognized\n".
	    "The options GFF2 or GFF3 are supported\n";
    }
}


#-----------------------------+
# PRINT REQUESTED HELP        |
#-----------------------------+
if ( ($show_usage) ) {
#    print_help ("usage", File::Spec->rel2abs($0) );
    print_help ("usage", $0 );
}

if ( ($show_help) || (!$ok) ) {
#    print_help ("help",  File::Spec->rel2abs($0) );
    print_help ("help",  $0 );
}

if ($show_man) {
    # User perldoc to generate the man documentation.
    system ("perldoc $0");
    exit($ok ? 0 : 2);
}

if ($show_version) {
    print "\nbatch_mask.pl:\n".
	"Version: $VERSION\n\n";
    exit;
}

my $in_gff = Bio::FeatureIO->new ('-file' => "$infile",
				  '-format' => 'GFF',
				  '-version' => 3);

my $x_in_gff = Bio::FeatureIO->new ('-file' => "$infile",
				  '-format' => 'GFF',
				  '-version' => 3);

#my $y_in_gff = Bio::FeatureIO->new ('-file' => "$infile",
#				  '-format' => 'GFF',
#				  '-version' => 3);

# unfortunately slow, and opening a new FILE for EACH
# 
my $overlap_count = 0;

# Object for graph, will initially do an undirected graph
# this is the scaffold graph
my $s_graph = Graph->new (directed => 0);

# First add vertices to the graph
my $gff_for_graph = Bio::FeatureIO->new ('-file' => "$infile",
					 '-format' => 'GFF',
					 '-version' => 3);

#-----------------------------------------------------------+
# WRITE NODES TO GRAPH                                      |
#-----------------------------------------------------------+
# Use GFF feature ID to tag the graph
# May need to add sanity check here to make ID usable by the 
# graph object
my $num_parent_features = 0;
while (my $g_feature = $gff_for_graph->next_feature() ) {

    if ($g_feature->primary_tag =~ $parent_feature) {  
	$s_graph->add_vertex( $g_feature->get_tag_values("ID") );
	$num_parent_features++;
	$all_sources{$g_feature->source_tag()} = 1;
#	print STDERR $g_feature->source_tag."\n";

    }
    
    
}


# Test of the gff overlaps method
while (my $x_feature = $x_in_gff->next_feature() ) {
    if ($x_feature->primary_tag =~ $parent_feature)  {
	print STDERR $x_feature->start."\n";
	print STDERR $x_feature->end."\n";
	# Get features that overlap the range of this feature
#	my @overlaping_features = $x_feature->features;
#	my @overlaping_features = $x_feature->overlaps;
    }
}

exit;

#-----------------------------------------------------------+
# WRITE EDGES TO GRAPH                                      |
#-----------------------------------------------------------+
while (my $x_feature = $x_in_gff->next_feature() ) {
    # Check that X is is primary tag
#    print STDERR $x_feature->primary_tag."\n";
    if ($x_feature->primary_tag =~ $parent_feature)  {

	my $y_in_gff = Bio::FeatureIO->new ('-file' => "$infile",
					    '-format' => 'GFF',
					    '-version' => 3);

	while (my $y_feature = $y_in_gff->next_feature() ) {

	    # Test for if both features are parents here
	    if ($y_feature->primary_tag =~ $parent_feature) {
		if ($x_feature->overlaps( $y_feature )) {
		    print STDERR "We have a overlap\n";
		    $overlap_count++;
		    if ($x_feature->equals($y_feature) ) {
			print STDERR "\tExact overlap\n";
		    }
		    my $int_range = $x_feature->intersection($y_feature);
		    my @x_ids = $x_feature->get_tag_values("ID");
		    # will assume a single ID
		    print STDERR "\tX_ID->".$x_ids[0]."\n";
		    my @y_ids = $y_feature->get_tag_values("ID");
		    print STDERR "\tY_ID".$y_ids[0]."\n";

		    # Add edge to graph object
		    my $x_node_label = $x_feature->get_tag_values("ID");
		    my $y_node_label = 	$y_feature->get_tag_values("ID");
		    $s_graph->add_edge( $x_node_label,
					$y_node_label);

		}
	    }
	} # End of iterating through y features
    }
    
}

#-----------------------------+
# DETERMINE CONNECTED         |
# COMPONENTS FOR GRAPH        |
#-----------------------------+
print STDERR "Determining the connected compoents graph\n";
my @cc = $s_graph->connected_components();

my $num_cc = @cc;

# The advantage of this appraoch is that can fetch subfeatures
# but this may require iterating through set to find ID match at a later as
# time
# @feats = $feat->sub_SeqFeature();

print STDERR "\n";
my @unique_sources = keys(%all_sources);
my $num_unique_sources = @unique_sources;

print STDERR "Number of unique sources: ".$num_unique_sources."\n";
print STDERR "Number of features (nodes): ".$num_parent_features."\n";
print STDERR "Number of overlaps (edges): ".$overlap_count."\n";
print STDERR "Number of Clusters: ".$num_cc."\n";


print STDERR "CLUSTER RESULTS:\n";
my $clust_num = 0;
foreach my $comp_list (@cc) {
    $clust_num++;
    my $num_clust_nodes = @$comp_list;
    for my $ind_comp (@$comp_list) {
	
    }
#    print 
}

exit;
__END__
# Early attempt at non-bioperl way below
#-----------------------------+
# OPEN FILEHANDLES            |
#-----------------------------+

if ($infile) {
    open (INFILE, "<$infile") ||
	die "Can not open input file\n";
}
else {
    open (INFILE, "<&STDIN") ||
	die "Can not open STDIN for input\n";
}


#-----------------------------+
# PROCESS GFF INPUT FILE      |
#-----------------------------+
my @parent_features;           # 2d array of parent features
my %all_sources = ();          # Sources array to extract unique list later
my $lin_num = 0;
while (<INFILE>) {
    chomp;
    $lin_num++;
    next if m/^#/;

    print STDERR "Processing $lin_num\n"
	if $verbose;
    
    my @gff_parts = split(/\t/, $_);
    my $num_parts = @gff_parts;
    unless ($num_parts == 9 ) {
	print STDERR "WARNING: Unexpected number of gff columns ".
	    "in line $lin_num. Skipping this line.\n"
    }

    
    if ($gff_parts[2] =~ $parent_feature) {
	print STDERR $_."\n"
	    if $verbose;

	push (@parent_features , \@gff_parts);
#	push (@all_sources, $gff_parts[1]);
	$all_sources{$gff_parts[1]} = 1;
	
    };

}

# Check that we have the parent features in an array
my $num_parent_features = @parent_features;
print STDERR "Num parent features $num_parent_features\n";

my @unique_sources = keys(%all_sources);
my $num_unique_sources = @unique_sources;
print STDERR "Num unique sources: ".$num_unique_sources."\n"
    if $verbose;
#-----------------------------+
# WORK ON GETTING LISTS FOR   |
# RANGE SET FUNCTIONS         |
#-----------------------------+


close (INFILE);


exit 0;

#-----------------------------------------------------------+ 
# SUBFUNCTIONS                                              |
#-----------------------------------------------------------+

sub print_help {
    my ($help_msg, $podfile) =  @_;
    # help_msg is the type of help msg to use (ie. help vs. usage)
    
    print "\n";
    
    #-----------------------------+
    # PIPE WITHIN PERL            |
    #-----------------------------+
    # This code made possible by:
    # http://www.perlmonks.org/index.pl?node_id=76409
    # Tie info developed on:
    # http://www.perlmonks.org/index.pl?node=perltie 
    #
    #my $podfile = $0;
    my $scalar = '';
    tie *STDOUT, 'IO::Scalar', \$scalar;
    
    if ($help_msg =~ "usage") {
	podselect({-sections => ["SYNOPSIS|MORE"]}, $0);
    }
    else {
	podselect({-sections => ["SYNOPSIS|ARGUMENTS|OPTIONS|MORE"]}, $0);
    }

    untie *STDOUT;
    # now $scalar contains the pod from $podfile you can see this below
    #print $scalar;

    my $pipe = IO::Pipe->new()
	or die "failed to create pipe: $!";
    
    my ($pid,$fd);

    if ( $pid = fork() ) { #parent
	open(TMPSTDIN, "<&STDIN")
	    or die "failed to dup stdin to tmp: $!";
	$pipe->reader();
	$fd = $pipe->fileno;
	open(STDIN, "<&=$fd")
	    or die "failed to dup \$fd to STDIN: $!";
	my $pod_txt = Pod::Text->new (sentence => 0, width => 78);
	$pod_txt->parse_from_filehandle;
	# END AT WORK HERE
	open(STDIN, "<&TMPSTDIN")
	    or die "failed to restore dup'ed stdin: $!";
    }
    else { #child
	$pipe->writer();
	$pipe->print($scalar);
	$pipe->close();	
	exit 0;
    }
    
    $pipe->close();
    close TMPSTDIN;

    print "\n";

    exit 0;
   
}

1;
__END__

=head1 NAME

Name.pl - Short program description. 

=head1 VERSION

This documentation refers to program version 0.1

=head1 SYNOPSIS

=head2 Usage

    Name.pl -i InFile -o OutFile

=head2 Required Arguments

    --infile        # Path to the input file
    --outfie        # Path to the output file

=head1 DESCRIPTION

This is what the program does

=head1 REQUIRED ARGUMENTS

=over 2

=item -i,--infile

Path of the input file.

=item -o,--outfile

Path of the output file.

=back

=head1 OPTIONS

=over 2

=item --usage

Short overview of how to use program from command line.

=item --help

Show program usage with summary of options.

=item --version

Show program version.

=item --man

Show the full program manual. This uses the perldoc command to print the 
POD documentation for the program.

=item -q,--quiet

Run the program with minimal output.

=back

=head1 EXAMPLES

The following are examples of how to use this script

=head2 Typical Use

This is a typcial use case.

=head1 DIAGNOSTICS

=over 2

=item * Expecting input from STDIN

If you see this message, it may indicate that you did not properly specify
the input sequence with -i or --infile flag. 

=back

=head1 CONFIGURATION AND ENVIRONMENT

Names and locations of config files
environmental variables
or properties that can be set.

=head1 DEPENDENCIES

=over

=item Data::Range::Compare

This perl module is required to find overlaps.

=back

Other modules or software that the program is dependent on.

=head1 BUGS AND LIMITATIONS

Any known bugs and limitations will be listed here.

=head1 REFERENCE

A manuscript is being submitted describing the DAWGPAWS program. 
Until this manuscript is published, please refer to the DAWGPAWS 
SourceForge website when describing your use of this program:

JC Estill and JL Bennetzen. 2009. 
The DAWGPAWS Pipeline for the Annotation of Genes and Transposable 
Elements in Plant Genomes.
http://dawgpaws.sourceforge.net/

=head1 LICENSE

GNU General Public License, Version 3

L<http://www.gnu.org/licenses/gpl.html>

=head1 AUTHOR

James C. Estill E<lt>JamesEstill at gmail.comE<gt>

=head1 HISTORY

STARTED:

UPDATED:

VERSION: $Rev$

=cut

#-----------------------------------------------------------+
# HISTORY                                                   |
#-----------------------------------------------------------+
#
