# -*-Perl-*- Test Harness script for Bioperl
# $Id: BioDBGFF.t 15407 2009-01-20 05:18:29Z cjfields $

use strict;
use Module::Build;

BEGIN {
    use lib '.';
    use Bio::Root::Test;
    
    test_begin(-tests => 279);
	
	use_ok('Bio::DB::GFF');
}

my $fasta_files = test_input_file('dbfa');
my $gff_file1   = test_input_file('biodbgff', 'test.gff');
my $gff_file2   = test_input_file('biodbgff', 'test.gff3');

my $build = Module::Build->current;
my $test_dsn = $build->notes('test_dsn');

my $adaptor = $test_dsn ? $test_dsn : 'memory';
$adaptor    = shift if @ARGV;

my @args;
if ($adaptor =~ /^dbi/) {
  my $cfg = {};
  $cfg->{dbd_driver} = $build->notes('dbd_driver');
  $cfg->{test_db} = $build->notes('test_db');
  $cfg->{test_host} = $build->notes('test_host');
  $cfg->{test_user} = $build->notes('test_user');
  $cfg->{test_pass} = $build->notes('test_pass');
  $cfg->{test_dsn} = $build->notes('test_dsn');
  
  $adaptor = "dbi::$cfg->{dbd_driver}" if $cfg->{dbd_driver};
  @args = ( '-adaptor'  => $adaptor,
	    '-dsn'     => $cfg->{test_dsn},
	  );
  push @args,('-user' => $cfg->{test_user}) if $cfg->{test_user};
  push @args,('-pass' => $cfg->{test_pass}) if $cfg->{test_pass};
} else {
  @args = ('-adaptor' => $adaptor,
	   '-create'  => 1);
}

push @args,('-aggregators' => ['transcript','processed_transcript']);

SKIP: {
for my $FILE ($gff_file1,$gff_file2) {

  my $db = eval { Bio::DB::GFF->new(@args) };
  skip "DB load failed? Skipping all! $@", 278 if $@;
  ok($db);

  $db->debug(0);
  $db->gff3_name_munging(1);

  # set the preferred groups
  $db->preferred_groups( [ 'transcript', 'gene', 'mRNA' ] );
  my @pg = $db->preferred_groups;
  is(scalar(@pg), 3);
  is($pg[1], 'gene'); 

  # exercise the loader
  ok($db->initialize(1));
  ok($db->load_gff($FILE));
  ok($db->load_fasta($fasta_files));

  # exercise db->types
  my @types = sort $db->types;
  is(scalar @types,11);
  is($types[0],'CDS:confirmed');
  is($types[-1],'transposon:tc1');
  my %types = $db->types('-enumerate'=>1);
  is($types{'transposon:tc1'},2);

  # exercise segment
  my $segment1 = $db->segment('Contig1');

  ok($segment1);
  is($segment1->length,37450);
  is($segment1->start,1);
  is($segment1->end,37450);
  is($segment1->strand,1);
  
  my $segment2  = $db->segment('Contig1',1=>1000);
  is($segment2->length,1000);
  is($segment2->start,1);
  is($segment2->end,1000);
  is($segment2->strand,1);
  
  my $segment3 = $db->segment('Contig1',10=>1);
  is($segment3->start,10);
  is($segment3->end,1);
  is($segment3->strand,-1);

  # exercise attribute fetching
  my @t = $db->fetch_feature_by_name(Transcript => 'trans-1');
  my ($t) = grep {$_->type eq 'transcript:confirmed'} @t;
  is($t->attributes('Note'),'function unknown');
  is(join(' ',sort $t->attributes('Gene')),'abc-1 xyz-2');
  my $att = $t->attributes;
  is(scalar @{$att->{Gene}},2);
  @t = sort {$a->display_name cmp $b->display_name} $db->fetch_feature_by_attribute('Gene'=>'abc-1');
  cmp_ok(@t,'>',0);
  is($t[0], $t);
  my $seg = $db->segment('Contig1');
  @t = $seg->features(-attributes=>{'Gene'=>'abc-1'});
  cmp_ok(@t,'>',0);
  is($seg->feature_count, 17);
  @t = $seg->features(-attributes=>{'Gene'=>'xyz-2',Note=>'Terribly interesting'});
  is(@t,1);

  # exercise dna() a bit
  my $dna = $segment2->dna;
  is(length $dna,1000);
  is(substr($dna,0,10),'gcctaagcct');
  is($segment3->dna,'aggcttaggc');
  is($segment1->dna, $db->dna($segment1->ref));

  # exercise ref()
  my $segment4 = $db->segment('-name'=>'c128.1','-class'=>'Transposon');
  is($segment4->length,1000);
  is($segment4->start,1);
  is($segment4->end,1000);
  is($segment4->ref,'c128.1');
  is($segment4->strand,1);
  ok(!$segment4->absolute);

  $segment4->absolute(1);
  ok($segment4->absolute);
  is($segment4->ref,'Contig1');
  is($segment4->start,5001);
  $segment4->absolute(0);
  my $tmp = $db->segment('Contig1',5001=>6000);
  is($segment4->dna,$tmp->dna);

  $segment4->ref('Contig1');
  is($segment4->ref,'Contig1');
  is($segment4->start,5001);
  is($segment4->end,6000);

  my $segment5 = $db->segment('-name'=>'c128.2','-class'=>'Transposon');
  is($segment5->length,1000);
  is($segment5->start,1);
  is($segment5->end,1000);
  is($segment5->ref,'c128.2');
  is($segment5->strand,1);

  $tmp = $db->segment('Contig1',9000,8001);
  is($segment5->dna,$tmp->dna);
  $segment5->absolute(1);
  is($segment5->strand,-1);

  # rel/rel addressing
  # first two positive strand features
  $segment4 = $db->segment('-name'=>'c128.1','-class'=>'Transposon');
  my $start4 = $segment4->abs_start;
  $segment5  = $db->segment('Transcript' => 'trans-1');
  my $start5 = $segment5->abs_start;
  $segment4->ref($segment5);
  is($segment4->strand,1);
  is($segment4->start,$start4-$start5+1);
  is($segment4->stop,$start4-$start5+$segment4->length);

  $segment4->ref('Transposon' => 'c128.1');
  $segment5->ref('Transcript' => 'trans-1');
  $segment5->ref($segment4);
  is($segment5->start,$start5-$start4+1);

  # now a positive on a negative strand feature
  my $segment6 = $db->segment('Transcript'=>'trans-2');
  my $start6 = $segment6->abs_start;
  is($segment6->strand,1);
  is($segment6->abs_strand,-1);
  $segment6->ref($segment4);
  is($segment6->start,$start6-$start4+1);
  is($segment6->strand,-1);

  $segment4->ref($segment6);
  is($segment4->start,$start6-$start4+1);
  is($segment4->strand,-1);
  is($segment4->ref,$segment6);

  # the reference sequence shouldn't affect the dna
  $segment6 = $db->segment('Transcript'=>'trans-2');
  $dna = $segment6->dna;
  $segment6->ref($segment4);
  is($segment6->dna,$dna);

  # segments should refuse to accept a reference sequence on a foreign segment
  undef $@;
  my $result = eval { $segment6->ref('Contig2') };
  ok(!$result);
  like($@, qr/are on different sequence segments/);

  # types across a segment
  $segment1 = $db->segment('Contig1');
  @types = sort $segment1->types;
  is(scalar @types,6);
  is($types[0],'CDS:confirmed');
  is($types[-1],'transposon:tc1');
  %types = $segment1->types('-enumerate'=>1);
  is($types{'similarity:est'},3);

  # features across a segment
  my @features = $segment1->features('-automerge'=>0);
  is(scalar @features,17);
  my %types_seen;
  foreach (@features) {
    $types_seen{$_->type}++;
  }
  my $inconsistency = 0;
  foreach (keys %types,keys %types_seen) {
    $inconsistency++ unless $types_seen{$_} == $types{$_};
  }
  ok(!$inconsistency);

  @features = sort {$a->start<=>$b->start} @features;
  is($features[0]->type,'Component:reference');
  is($features[-1]->type,'exon:confirmed');

  # make sure that we can use features to get at dna
  is($features[0]->dna,$db->segment('Contig1',$features[0]->start,$features[0]->end)->dna);

  # check three forward features and three reverse features
  # (This depends on the test.gff data)
  for (1..3,-3..-1) {
    $segment2 = $db->segment($features[$_],50,100);
    if ($features[$_]->strand >= 0) {
      is($segment2->dna,$db->segment('Contig1',
				     $features[$_]->start+50-1,
				     $features[$_]->start+100-1)->dna)
    } else {
      is($segment2->dna,$db->segment('Contig1',
				     $features[$_]->start-50+1,
				     $features[$_]->start-100+1)->dna)
    }
  }

  # exercise the aggregator
  my $aggregator = Bio::DB::GFF::Aggregator->new('-method'      => 'aggregated_transcript',
						 '-main_method' => 'transcript',
						 '-sub_parts'   => ['exon','CDS']);
  $db->add_aggregator($aggregator);
  $segment1 = $db->segment('Contig1');
  @features = sort $segment1->features('aggregated_transcript');   # sort so that trans-1 comes first
  is(scalar @features,2);
  cmp_ok($features[0]->Exon, '>', 0);
  cmp_ok($features[0]->Cds,'>', 0);

  # Test that sorting is correct.  The way that test.gff is set up, the lower one is
  # on the + strand and the higher is on the -.
  @features = sort {$a->start <=> $b->start} @features;
  is($features[0]->strand,1);
  is($features[1]->strand,-1);

  my $last = 0;
  $inconsistency = 0;
  foreach ($features[0]->Exon) {
    $inconsistency++ if $_->start > $_->end;
    $inconsistency++ if $last && $_->start < $last;
    $last = $_->start;
  }
  ok(!$inconsistency);

  $inconsistency = $last = 0;
  foreach ($features[1]->Exon) {
    $inconsistency++ if $_->start < $_->end;
    $inconsistency++ if $last && $_->start > $last;
    $last = $_->start;
  }
  ok(!$inconsistency);
  
  # relative addressing in aggregated features
  my $transcript1 = $db->segment($features[0]);
  $transcript1->ref($features[0]);
  my @overlap     = sort {$a->start <=> $b->start } $transcript1->features;
  is(scalar(@overlap),5);
  is($overlap[0]->start,-999);

  $transcript1 = $db->segment('Transcript' => 'trans-1');
  @overlap     = sort {$a->start <=> $b->start } $transcript1->features;
  is($overlap[0]->start,-999);

  # test strandedness of features
  $segment1 = $db->segment('-class' => 'Transcript',
			   '-name'  => 'trans-3',
			   '-start' => 1,
			   '-stop'  => 6000);
  is($segment1->strand,1);
  @overlap  = sort {$a->start <=> $b->start} $segment1->features('transcript');
  is(scalar(@overlap),2);
  is($overlap[0]->name,'trans-3');
  is($overlap[1]->name,'trans-4');
  is($overlap[0]->strand,1);
  is($overlap[1]->strand,-1);

  # testing feature id and group_id
  my $tf = $overlap[0];
  ok(defined $tf->id);
  my $t1 = $db->fetch_feature_by_id($tf->id);
  is($t1->id,$tf->id);

  SKIP: {
    if (defined $tf->group_id) {
      my $t2 = $db->fetch_feature_by_gid($tf->group_id);
      is($t2->group_id,$tf->group_id);
      is($t2->group_id,$t1->group_id);
    } else {
      skip("fetch_feature_by_gid() not implemented by this adaptor",2);
    }
  }

  $segment1 = $db->segment('-class' => 'Transcript',
			   '-name'  => 'trans-4',
			   '-start' => 1,
			   '-stop'  => 6000);
  is($segment1->strand,1);
  @overlap = sort {$a->start <=> $b->start} $segment1->features('transcript');
  is($overlap[0]->name,'trans-4');
  is($overlap[1]->name,'trans-3');
  is($overlap[0]->strand,1);
  is($overlap[1]->strand,-1);

  @overlap = sort {$a->start <=> $b->start} $segment1->features('Component');
  is($overlap[0]->strand,0);

SKIP: {
  # test preferred group assignments
  if ($FILE =~ /\.gff$/) {
    my @gene = $db->get_feature_by_name( gene => 'gene-9' );
    my @mrna = $db->get_feature_by_name( mRNA => 'trans-9' );
    is($gene[0]->ref, 'Contig4');
    is(scalar(@gene), 2);
    is(scalar(@mrna), 1);
  } else {
    skip('preferred groups are not supported by gff3',3);
  }
}

  # test iterator across a segment
  $segment1 = $db->segment('Contig1');
  my $i = $segment1->features('-automerge'=>0,'-iterator'=>1);
  my %strand;
  while (my $s = $i->next_feature) {
    $strand{$s->strand}++;
  }
  is(keys %strand, 3);

  # test iterator across entire database
  $i = $db->features('-automerge'=>0,'-iterator'=>1);
  %strand = ();
  while (my $s = $i->next_feature) {
    $strand{$s->strand}++;
  }
  is(keys %strand, 3);

  # test iterator across a segment, limited by an attribute
  $i = $seg->get_feature_stream(-attributes=>{'Gene'=>'abc-1',Note=>'function unknown'});
  my $count = 0;
  while ($i->next_seq) {
    $count++;
  }
  is($count,2);

  # test that aliases work
  my $st1 = $db->segment(Transcript => 'trans-3');
  ok($st1);
  my $st2 = $db->segment(Transcript => 'trans-18');  # this is an alias!
  ok($st2);
  is($st1,$st2);
  my @transcripts = $st1->features('transcript');
  is(($transcripts[0]->aliases)[0],'trans-18');

  # test truncation
  $db->strict_bounds_checking(1);
  my $tseg = $db->segment(-name=>'trans-1',-class=>'Transcript',-start=>1,-stop=>500);
  ok(!$tseg->truncated);
  $tseg    = $db->segment(-name=>'trans-1',-class=>'Transcript',-start=>1,-stop=>50000);
  ok($tseg->truncated);
  $db->strict_bounds_checking(0);
  $tseg    = $db->segment(-name=>'trans-1',-class=>'Transcript',-start=>1,-stop=>50000);
  ok(!$tseg->truncated);

  # test the processed_transcript aggregator
  $db->clear_aggregators;
  $db->add_aggregator('processed_transcript');
  my @f = $db->fetch_feature_by_name(mRNA => 'trans-8');
  is(scalar @f,1);
  is($f[0]->length,35000-32000+1);
  is(scalar $f[0]->CDS,3);
  is(scalar $f[0]->UTR,2);

  # test deletions
  # segment delete() method
  my $clone = $db->segment(Clone=>'M7.3');
  my $overlapping_feature_count = $clone->features(-range_type =>'overlaps');
  my $contained_feature_count   = $clone->features(-range_type =>'contains');
  is(scalar $clone->delete(-range_type=>'contains'),$contained_feature_count);
  is(scalar $clone->features,$overlapping_feature_count - $contained_feature_count);

  # database delete() method
  is($db->delete(-type=>['mRNA:confirmed','transposon:tc1']),4);
  is($db->delete(-type=>'UTR',-ref=>'Contig29'),undef);
  is($db->delete(-type=>'CDS',-ref=>'AL12345.2',-class=>'Clone'),3);
  is($db->delete_features(1,2,3),3);

  SKIP: {
    $result = eval {
      is($db->delete_groups(1,2,3,4,5),5);
      my @features = $db->get_feature_by_name(Sequence => 'Contig2');
      is($db->delete_groups(@features),1);
      1;
    };
    if (!$result && $@ =~ /not implemented/i) {
      skip("delete_groups() not implemented by this adaptor",2);
    }
  }
  
  SKIP: {
	test_skip(-tests => 1, -excludes_os => 'mswin');
	
	# test ability to pass adaptors across a fork
	if (my $child = open(F,"-|")) { # parent reads from child
		ok(scalar <F>);
		close F;
	}
	else { # in child
		$db->clone;
		my @f = $db->features();
		print @f>0;
		exit 0;
	}
  }

  ok(!defined eval{$db->delete()});
  ok($db->delete(-force=>1));
  is(scalar $db->features,0);
  ok(!$db->segment('Contig1'));

}

}



END {
  unlink $fasta_files."/directory.index";
}