#!/usr/local/bin/perl -w
use strict;

use Bio::Chaos::ChaosGraph;
use Bio::SeqIO;
use Getopt::Long;
use FileHandle;

my ($fmt, $outfmt, $type, $writer) =
  qw(genbank chaos seq xml);
my @remove_types = ();
my $make_islands;
my $ascii;
my $nameby = 'feature_id';
my $error_threshold = 3;
my $dir = '.';
my $ds_root;
my $ds;    # The pathname to the datastore root.
my $verbose;  # reference to the datastore object.
my $include_haplotypes;
GetOptions("fmt|i=s"=>\$fmt,
           "writer|w=s"=>\$writer,
           "outfmt|o=s"=>\$outfmt,
           "type|t=s"=>\$type,
	   "islands"=>\$make_islands,
           "dir=s"=>\$dir,
	   "ascii|a"=>\$ascii,
	   "nameby|n=s"=>\$nameby,
	   "ethresh|e=s"=>\$error_threshold,
	   "remove_type=s@"=>\@remove_types,
           "ds_root=s"=>\$ds_root,
           "verbose|v=s"=>\$verbose,
           "include_haplotypes"=>\$include_haplotypes,
	   "help|h"=>sub {
	       system("perldoc $0"); exit 0;
	   }
	  );

my $W = Data::Stag->getformathandler($writer);
$W->fh(\*STDOUT);

if ($ds_root) {
    require "Datastore/MD5.pm";
    $ds = new Datastore::MD5("root" => $ds_root, "depth" => 2);
}
else {
    mkdir($dir) unless $dir eq '.';
}

foreach my $file (@ARGV) {
    debug("converting $file") if $verbose;
    my $chaos;
    my $infh;
    if ($file =~ /\.gz$/) {
        $infh = FileHandle->new("gzip -dc $file |");
    }
    else {
        $infh = FileHandle->new($file);
    }
    if (!$infh) {
        print STDERR "Cannot open: $file for reading\n";
        exit 1;
    }
    my $seqio =
      Bio::SeqIO->new(-fh=> $infh,
                      -format => $fmt);
    while (my $seq = $seqio->next_seq()) {
        if ($seq->desc =~ /haplotype/i) {
            # this is NASTY, but unfortunately GB has no
            # consistent way of defining haps;
            # in general we wish to remove alternative haps 
            # an only show the reference sequence
            next unless $include_haplotypes;
        }
        debug("Got bioperl Seq: %s %d bp [$file]\n", $seq->accession, $seq->length) if $verbose;
        my $C =
          Bio::Chaos::ChaosGraph->new;
        $C->verbose($verbose);
        my $unflattener = $C->unflattener;
        my $type_mapper = $C->type_mapper;
        $unflattener->error_threshold($error_threshold);
        $unflattener->verbose($verbose) if $verbose;
        if (@remove_types) {
            debug("removing types %s", join(' ',@remove_types)) if $verbose;
            $unflattener->remove_types(-seq=>$seq,
                                       -types=>\@remove_types);
        }
        debug("Unflattening Seq: %s [$file]\n", $seq->accession) if $verbose;
        eval {
            $unflattener->unflatten_seq(-seq=>$seq,
                                        -use_magic=>1);
            $unflattener->report_problems(\*STDERR);
            debug("Mapping types") if $verbose;
            $type_mapper->map_types_to_SO(-seq=>$seq);
        };
        if ($@) {
            print $@;
            printf STDERR "  Problems unflattening: %s BYE BYE [$file]\n", $seq->accession;
            exit 1;
        }
        debug("Unflattened Seq: %s [$file]\n", $seq->accession) if $verbose;

        my $outio = Bio::SeqIO->new( -format => 'chaos');
        $outio->write_seq($seq); # "writes" to a stag object
        $outio->end_of_data;

        # free memory
        %$seq = ();
        my $stag = $outio->handler->stag;
        debug("initiating chaos graph") if $verbose;
        $C->init_from_stag($stag);
	
        my $seqfs = $C->top_unlocalised_features;
        if (@$seqfs == 0) {
            $C->freak("no top level feature");
        }
        if (@$seqfs > 1) {
            $C->freak("top unlocalised feats!=1", @$seqfs);
        }
        my $seqf = shift @$seqfs;
        my $islands;
        if ($make_islands) {
            debug("making islands") if $verbose;
            my $acc = $seqf->get_feature_id;
            my $path_prefix = "$dir/$acc";
            mkdir($path_prefix) unless -d $path_prefix;

            my $fs = $C->top_features;
            my @islands = ();
            foreach my $f (@$fs) {
                my $type = $f->get_type;
                $C->freak("no type", $f) unless $type;
                next unless $f->get_type eq 'gene';
                eval {
                    debug("Making island around: %s\n", $f->get_feature_id)
                      if $verbose;
                    my $islandC = $C->make_island($f, 500);
                    #		print $islandC->asciitree;
                    my $W = Data::Stag->getformathandler($writer);
                    my $id = $f->get($nameby);
                    $id =~ tr/A-Za-z0-9_:;\.//cd; # make name safe
                    my $fh;
                    my $fn;
                    if ($ds_root) {
                        $ds->mkdir($id) || die "Unable to ds->mkdir() for $id";
                        $fn = sprintf("%s/%s.%s",
                                      $ds->id_to_dir($id),
                                      $id,
                                      $writer);
                    }
                    else {
                        $fn = sprintf("%s/%s.%s",
                                      $path_prefix,
                                      $id,
                                      $writer);
                    }
                    $fh = FileHandle->new(">$fn") || die "can't open $fn";
                    $W->fh($fh);

                    # temporarily clear seq residues
                    # (the island contig will have the local residues)
                    my $res = $seqf->sget_residues;
                    $seqf->unset_residues;
                    $islandC->metadata(Data::Stag->new(focus_feature_id=>$f->get_feature_id));
                    debug("writing gene island file %s", $fn)
                      if $verbose;
                    $islandC->stag->events($W);
                    $seqf->set_residues($res);
                    $fh->close;
                };
                if ($@) {
                    print STDERR $@;
                    print STDERR "Problem resides with source feature\n";
                    print STDERR $f->xml;
                    exit 1;
                }
            }
        } else {
            #	print $C->asciitree;
            if ($ascii) {
                print $C->asciitree;
            } else {
                $chaos = $C->stag;
                $chaos->sax($W);
            }
        }

    }
    debug("$file converted") if $verbose;
    $infh->close;
}
debug("ALL DONE!") if $verbose;

exit 0;

# subs
sub debug {
    return unless $verbose;
    my $fmt = shift;
    my $t = time;
    my $lt = localtime $t;
    print STDERR "# "; 
    printf STDERR ($fmt, @_);
    print STDERR " [$t] $lt\n";
}

__END__

=head1 NAME 

  cx-genbank2chaos.pl.pl

=head1 SYNOPSIS

  cx-genbank2chaos.pl.pl sample-data/AE003734.gbk > AE003734.chaos.xml

  cx-genbank2chaos.pl.pl -islands sample-data/AE003734.gbk

=head1 DESCRIPTION

Converts a genbank file to a chaos xml file (or a collection of chaos
xml files).

The genbank file is 'unflattened' in order to infer the relationships
between features

with the -islands option set, this loops through a list of
genbank-formatted files and builds a chaos file for every gene

by default it will store each gene in a directory named by the
sequence accession. it will name each file by the unique feature_id;
for example

  AE003644.2/
    gene:EMBLGenBankSwissProt:AE003644:128108:128179.xml
    gene:EMBLGenBankSwissProt:AE003644:128645:128716.xml
    gene:EMBLGenBankSwissProt:AE003644:128923:128994.xml

You can change the field used to name the file with -nameby; for
example, if you use the chado/chaos B<name> field like this:

  cx-genbank2chaos.pl.pl -islands -nameby name AE003734.gbk

You will get

  AE003644.3/
   noc.xml
   osp.xml
   BG:DS07721.3.xml

the default is the B<feature_id> field, which is usually more
unix-friendly (fly genes can have all kinds of weird characters in
their name); also using the 'name' field could run into uniqueness
issues.

=head1 HOW IT WORKS

=over

=item 1 - parse genbank to bioperl

uses L<Bio::SeqIO::genbank>

=item 2 - unflatten the flat list of bioperl SeqFeatures

uses L<Bio::Seqfeature::Tools::Unflattener>

=item 3 - turn bioperl objects into chaos datastructure

uses L<Bio::SeqIO::chaos>

=item 4 - remap every gene to an 'island' (virtual contig)

uses L<Bio::Chaos::ChaosGraph>

=item 5 - spit out each virtual contig chaos graph to a file

uses L<Bio::Chaos::ChaosGraph>

=back

=head1 ARGUMENTS

=head2 -islands

exports one file per gene

=head2 -ethresh ERRORTHRESH

Sets the error threshold. See L<Bio::SeqFeature::Tools::Unflattener>

you will want to keep this at its default setting of 3 (insensitive)

=head2 -remove_type GENBANKFEATURETYPE

This will remove all features of a certain type prior to unflattening

This is useful if you wish to exclude a certain kind of feature (eg
variation) from your analysis

It is also required for the genbank release of S_Pombe, which has a
few scattered types purportedly of mRNA which confuse the unflattening
process

=head2 -ds_root DIR b<EXPERIMENTAL>

Root directory for building a datastore - see L<Datastore::MD5>

=head2 -include_haplotypes

by default, only reference sequences are exported. if the genbank
definition like contains the string "haplotype", then this is probably
an alternative haplotype that will skew analyses. this is removed by
default, unless this switch is set

For an example, see contigs NG_002432 and NT_007592 (the former is an
alternate haplotype of the latter)


=head1 REQUIREMENTS

bioperl 1.5 or later

=cut