t/LocalDB/BioDBGFF.t

   1 # -*-Perl-*- Test Harness script for Bioperl
   2 # $Id$
   3
   4 use strict;
   5 use Module::Build;
   6
   7 BEGIN {
   8     use lib '.';
   9     use Bio::Root::Test;
  10
  11     test_begin(-tests => 279);
  12
  13         use_ok('Bio::DB::GFF');
  14 }
  15
  16 my $fasta_files = test_input_file('dbfa');
  17 my $gff_file1   = test_input_file('biodbgff', 'test.gff');
  18 my $gff_file2   = test_input_file('biodbgff', 'test.gff3');
  19
  20 my $build = Module::Build->current;
  21 my $test_dsn = $build->notes('test_dsn');
  22
  23 my $adaptor = $test_dsn ? $test_dsn : 'memory';
  24 $adaptor    = shift if @ARGV;
  25
  26 my @args;
  27 if ($adaptor =~ /^dbi/) {
  28   my $cfg = {};
  29   $cfg->{dbd_driver} = $build->notes('dbd_driver');
  30   $cfg->{test_db} = $build->notes('test_db');
  31   $cfg->{test_host} = $build->notes('test_host');
  32   $cfg->{test_user} = $build->notes('test_user');
  33   $cfg->{test_pass} = $build->notes('test_pass');
  34   $cfg->{test_dsn} = $build->notes('test_dsn');
  35
  36   $adaptor = "dbi::$cfg->{dbd_driver}" if $cfg->{dbd_driver};
  37   @args = ( '-adaptor'  => $adaptor,
  38             '-dsn'     => $cfg->{test_dsn},
  39           );
  40   push @args,('-user' => $cfg->{test_user}) if $cfg->{test_user};
  41   push @args,('-pass' => $cfg->{test_pass}) if $cfg->{test_pass};
  42 } else {
  43   @args = ('-adaptor' => $adaptor,
  44            '-create'  => 1);
  45 }
  46
  47 push @args,('-aggregators' => ['transcript','processed_transcript']);
  48
  49 SKIP: {
  50 for my $FILE ($gff_file1,$gff_file2) {
  51
  52   my $db = eval { Bio::DB::GFF->new(@args) };
  53   skip "DB load failed? Skipping all! $@", 278 if $@;
  54   ok($db);
  55
  56   $db->debug(0);
  57   $db->gff3_name_munging(1);
  58
  59   # set the preferred groups
  60   $db->preferred_groups( [ 'transcript', 'gene', 'mRNA' ] );
  61   my @pg = $db->preferred_groups;
  62   is(scalar(@pg), 3);
  63   is($pg[1], 'gene');
  64
  65   # exercise the loader
  66   ok($db->initialize(1));
  67   ok($db->load_gff($FILE));
  68   ok($db->load_fasta($fasta_files));
  69
  70   # exercise db->types
  71   my @types = sort $db->types;
  72   is(scalar @types,11);
  73   is($types[0],'CDS:confirmed');
  74   is($types[-1],'transposon:tc1');
  75   my %types = $db->types('-enumerate'=>1);
  76   is($types{'transposon:tc1'},2);
  77
  78   # exercise segment
  79   my $segment1 = $db->segment('Contig1');
  80
  81   ok($segment1);
  82   is($segment1->length,37450);
  83   is($segment1->start,1);
  84   is($segment1->end,37450);
  85   is($segment1->strand,1);
  86
  87   my $segment2  = $db->segment('Contig1',1=>1000);
  88   is($segment2->length,1000);
  89   is($segment2->start,1);
  90   is($segment2->end,1000);
  91   is($segment2->strand,1);
  92
  93   my $segment3 = $db->segment('Contig1',10=>1);
  94   is($segment3->start,10);
  95   is($segment3->end,1);
  96   is($segment3->strand,-1);
  97
  98   # exercise attribute fetching
  99   my @t = $db->fetch_feature_by_name(Transcript => 'trans-1');
 100   my ($t) = grep {$_->type eq 'transcript:confirmed'} @t;
 101   is($t->attributes('Note'),'function unknown');
 102   is(join(' ',sort $t->attributes('Gene')),'abc-1 xyz-2');
 103   my $att = $t->attributes;
 104   is(scalar @{$att->{Gene}},2);
 105   @t = sort {$a->display_name cmp $b->display_name} $db->fetch_feature_by_attribute('Gene'=>'abc-1');
 106   cmp_ok(@t,'>',0);
 107   is($t[0], $t);
 108   my $seg = $db->segment('Contig1');
 109   @t = $seg->features(-attributes=>{'Gene'=>'abc-1'});
 110   cmp_ok(@t,'>',0);
 111   is($seg->feature_count, 17);
 112   @t = $seg->features(-attributes=>{'Gene'=>'xyz-2',Note=>'Terribly interesting'});
 113   is(@t,1);
 114
 115   # exercise dna() a bit
 116   my $dna = $segment2->dna;
 117   is(length $dna,1000);
 118   is(substr($dna,0,10),'gcctaagcct');
 119   is($segment3->dna,'aggcttaggc');
 120   is($segment1->dna, $db->dna($segment1->ref));
 121
 122   # exercise ref()
 123   my $segment4 = $db->segment('-name'=>'c128.1','-class'=>'Transposon');
 124   is($segment4->length,1000);
 125   is($segment4->start,1);
 126   is($segment4->end,1000);
 127   is($segment4->ref,'c128.1');
 128   is($segment4->strand,1);
 129   ok(!$segment4->absolute);
 130
 131   $segment4->absolute(1);
 132   ok($segment4->absolute);
 133   is($segment4->ref,'Contig1');
 134   is($segment4->start,5001);
 135   $segment4->absolute(0);
 136   my $tmp = $db->segment('Contig1',5001=>6000);
 137   is($segment4->dna,$tmp->dna);
 138
 139   $segment4->ref('Contig1');
 140   is($segment4->ref,'Contig1');
 141   is($segment4->start,5001);
 142   is($segment4->end,6000);
 143
 144   my $segment5 = $db->segment('-name'=>'c128.2','-class'=>'Transposon');
 145   is($segment5->length,1000);
 146   is($segment5->start,1);
 147   is($segment5->end,1000);
 148   is($segment5->ref,'c128.2');
 149   is($segment5->strand,1);
 150
 151   $tmp = $db->segment('Contig1',9000,8001);
 152   is($segment5->dna,$tmp->dna);
 153   $segment5->absolute(1);
 154   is($segment5->strand,-1);
 155
 156   # rel/rel addressing
 157   # first two positive strand features
 158   $segment4 = $db->segment('-name'=>'c128.1','-class'=>'Transposon');
 159   my $start4 = $segment4->abs_start;
 160   $segment5  = $db->segment('Transcript' => 'trans-1');
 161   my $start5 = $segment5->abs_start;
 162   $segment4->ref($segment5);
 163   is($segment4->strand,1);
 164   is($segment4->start,$start4-$start5+1);
 165   is($segment4->stop,$start4-$start5+$segment4->length);
 166
 167   $segment4->ref('Transposon' => 'c128.1');
 168   $segment5->ref('Transcript' => 'trans-1');
 169   $segment5->ref($segment4);
 170   is($segment5->start,$start5-$start4+1);
 171
 172   # now a positive on a negative strand feature
 173   my $segment6 = $db->segment('Transcript'=>'trans-2');
 174   my $start6 = $segment6->abs_start;
 175   is($segment6->strand,1);
 176   is($segment6->abs_strand,-1);
 177   $segment6->ref($segment4);
 178   is($segment6->start,$start6-$start4+1);
 179   is($segment6->strand,-1);
 180
 181   $segment4->ref($segment6);
 182   is($segment4->start,$start6-$start4+1);
 183   is($segment4->strand,-1);
 184   is($segment4->ref,$segment6);
 185
 186   # the reference sequence shouldn't affect the dna
 187   $segment6 = $db->segment('Transcript'=>'trans-2');
 188   $dna = $segment6->dna;
 189   $segment6->ref($segment4);
 190   is($segment6->dna,$dna);
 191
 192   # segments should refuse to accept a reference sequence on a foreign segment
 193   undef $@;
 194   my $result = eval { $segment6->ref('Contig2') };
 195   ok(!$result);
 196   like($@, qr/are on different sequence segments/);
 197
 198   # types across a segment
 199   $segment1 = $db->segment('Contig1');
 200   @types = sort $segment1->types;
 201   is(scalar @types,6);
 202   is($types[0],'CDS:confirmed');
 203   is($types[-1],'transposon:tc1');
 204   %types = $segment1->types('-enumerate'=>1);
 205   is($types{'similarity:est'},3);
 206
 207   # features across a segment
 208   my @features = $segment1->features('-automerge'=>0);
 209   is(scalar @features,17);
 210   my %types_seen;
 211   foreach (@features) {
 212     $types_seen{$_->type}++;
 213   }
 214   my $inconsistency = 0;
 215   foreach (keys %types,keys %types_seen) {
 216     $inconsistency++ unless $types_seen{$_} == $types{$_};
 217   }
 218   ok(!$inconsistency);
 219
 220   @features = sort {$a->start<=>$b->start} @features;
 221   is($features[0]->type,'Component:reference');
 222   is($features[-1]->type,'exon:confirmed');
 223
 224   # make sure that we can use features to get at dna
 225   is($features[0]->dna,$db->segment('Contig1',$features[0]->start,$features[0]->end)->dna);
 226
 227   # check three forward features and three reverse features
 228   # (This depends on the test.gff data)
 229   for (1..3,-3..-1) {
 230     $segment2 = $db->segment($features[$_],50,100);
 231     if ($features[$_]->strand >= 0) {
 232       is($segment2->dna,$db->segment('Contig1',
 233                                      $features[$_]->start+50-1,
 234                                      $features[$_]->start+100-1)->dna)
 235     } else {
 236       is($segment2->dna,$db->segment('Contig1',
 237                                      $features[$_]->start-50+1,
 238                                      $features[$_]->start-100+1)->dna)
 239     }
 240   }
 241
 242   # exercise the aggregator
 243   my $aggregator = Bio::DB::GFF::Aggregator->new('-method'      => 'aggregated_transcript',
 244                                                  '-main_method' => 'transcript',
 245                                                  '-sub_parts'   => ['exon','CDS']);
 246   $db->add_aggregator($aggregator);
 247   $segment1 = $db->segment('Contig1');
 248   @features = sort $segment1->features('aggregated_transcript');   # sort so that trans-1 comes first
 249   is(scalar @features,2);
 250   cmp_ok($features[0]->Exon, '>', 0);
 251   cmp_ok($features[0]->Cds,'>', 0);
 252
 253   # Test that sorting is correct.  The way that test.gff is set up, the lower one is
 254   # on the + strand and the higher is on the -.
 255   @features = sort {$a->start <=> $b->start} @features;
 256   is($features[0]->strand,1);
 257   is($features[1]->strand,-1);
 258
 259   my $last = 0;
 260   $inconsistency = 0;
 261   foreach ($features[0]->Exon) {
 262     $inconsistency++ if $_->start > $_->end;
 263     $inconsistency++ if $last && $_->start < $last;
 264     $last = $_->start;
 265   }
 266   ok(!$inconsistency);
 267
 268   $inconsistency = $last = 0;
 269   foreach ($features[1]->Exon) {
 270     $inconsistency++ if $_->start < $_->end;
 271     $inconsistency++ if $last && $_->start > $last;
 272     $last = $_->start;
 273   }
 274   ok(!$inconsistency);
 275
 276   # relative addressing in aggregated features
 277   my $transcript1 = $db->segment($features[0]);
 278   $transcript1->ref($features[0]);
 279   my @overlap     = sort {$a->start <=> $b->start } $transcript1->features;
 280   is(scalar(@overlap),5);
 281   is($overlap[0]->start,-999);
 282
 283   $transcript1 = $db->segment('Transcript' => 'trans-1');
 284   @overlap     = sort {$a->start <=> $b->start } $transcript1->features;
 285   is($overlap[0]->start,-999);
 286
 287   # test strandedness of features
 288   $segment1 = $db->segment('-class' => 'Transcript',
 289                            '-name'  => 'trans-3',
 290                            '-start' => 1,
 291                            '-stop'  => 6000);
 292   is($segment1->strand,1);
 293   @overlap  = sort {$a->start <=> $b->start} $segment1->features('transcript');
 294   is(scalar(@overlap),2);
 295   is($overlap[0]->name,'trans-3');
 296   is($overlap[1]->name,'trans-4');
 297   is($overlap[0]->strand,1);
 298   is($overlap[1]->strand,-1);
 299
 300   # testing feature id and group_id
 301   my $tf = $overlap[0];
 302   ok(defined $tf->id);
 303   my $t1 = $db->fetch_feature_by_id($tf->id);
 304   is($t1->id,$tf->id);
 305
 306   SKIP: {
 307     if (defined $tf->group_id) {
 308       my $t2 = $db->fetch_feature_by_gid($tf->group_id);
 309       is($t2->group_id,$tf->group_id);
 310       is($t2->group_id,$t1->group_id);
 311     } else {
 312       skip("fetch_feature_by_gid() not implemented by this adaptor",2);
 313     }
 314   }
 315
 316   $segment1 = $db->segment('-class' => 'Transcript',
 317                            '-name'  => 'trans-4',
 318                            '-start' => 1,
 319                            '-stop'  => 6000);
 320   is($segment1->strand,1);
 321   @overlap = sort {$a->start <=> $b->start} $segment1->features('transcript');
 322   is($overlap[0]->name,'trans-4');
 323   is($overlap[1]->name,'trans-3');
 324   is($overlap[0]->strand,1);
 325   is($overlap[1]->strand,-1);
 326
 327   @overlap = sort {$a->start <=> $b->start} $segment1->features('Component');
 328   is($overlap[0]->strand,0);
 329
 330 SKIP: {
 331   # test preferred group assignments
 332   if ($FILE =~ /\.gff$/) {
 333     my @gene = $db->get_feature_by_name( gene => 'gene-9' );
 334     my @mrna = $db->get_feature_by_name( mRNA => 'trans-9' );
 335     is($gene[0]->ref, 'Contig4');
 336     is(scalar(@gene), 2);
 337     is(scalar(@mrna), 1);
 338   } else {
 339     skip('preferred groups are not supported by gff3',3);
 340   }
 341 }
 342
 343   # test iterator across a segment
 344   $segment1 = $db->segment('Contig1');
 345   my $i = $segment1->features('-automerge'=>0,'-iterator'=>1);
 346   my %strand;
 347   while (my $s = $i->next_feature) {
 348     $strand{$s->strand}++;
 349   }
 350   is(keys %strand, 3);
 351
 352   # test iterator across entire database
 353   $i = $db->features('-automerge'=>0,'-iterator'=>1);
 354   %strand = ();
 355   while (my $s = $i->next_feature) {
 356     $strand{$s->strand}++;
 357   }
 358   is(keys %strand, 3);
 359
 360   # test iterator across a segment, limited by an attribute
 361   $i = $seg->get_feature_stream(-attributes=>{'Gene'=>'abc-1',Note=>'function unknown'});
 362   my $count = 0;
 363   while ($i->next_seq) {
 364     $count++;
 365   }
 366   is($count,2);
 367
 368   # test that aliases work
 369   my $st1 = $db->segment(Transcript => 'trans-3');
 370   ok($st1);
 371   my $st2 = $db->segment(Transcript => 'trans-18');  # this is an alias!
 372   ok($st2);
 373   is($st1,$st2);
 374   my @transcripts = $st1->features('transcript');
 375   is(($transcripts[0]->aliases)[0],'trans-18');
 376
 377   # test truncation
 378   $db->strict_bounds_checking(1);
 379   my $tseg = $db->segment(-name=>'trans-1',-class=>'Transcript',-start=>1,-stop=>500);
 380   ok(!$tseg->truncated);
 381   $tseg    = $db->segment(-name=>'trans-1',-class=>'Transcript',-start=>1,-stop=>50000);
 382   ok($tseg->truncated);
 383   $db->strict_bounds_checking(0);
 384   $tseg    = $db->segment(-name=>'trans-1',-class=>'Transcript',-start=>1,-stop=>50000);
 385   ok(!$tseg->truncated);
 386
 387   # test the processed_transcript aggregator
 388   $db->clear_aggregators;
 389   $db->add_aggregator('processed_transcript');
 390   my @f = $db->fetch_feature_by_name(mRNA => 'trans-8');
 391   is(scalar @f,1);
 392   is($f[0]->length,35000-32000+1);
 393   is(scalar $f[0]->CDS,3);
 394   is(scalar $f[0]->UTR,2);
 395
 396   # test deletions
 397   # segment delete() method
 398   my $clone = $db->segment(Clone=>'M7.3');
 399   my $overlapping_feature_count = $clone->features(-range_type =>'overlaps');
 400   my $contained_feature_count   = $clone->features(-range_type =>'contains');
 401   is(scalar $clone->delete(-range_type=>'contains'),$contained_feature_count);
 402   is(scalar $clone->features,$overlapping_feature_count - $contained_feature_count);
 403
 404   # database delete() method
 405   is($db->delete(-type=>['mRNA:confirmed','transposon:tc1']),4);
 406   is($db->delete(-type=>'UTR',-ref=>'Contig29'),undef);
 407   is($db->delete(-type=>'CDS',-ref=>'AL12345.2',-class=>'Clone'),3);
 408   is($db->delete_features(1,2,3),3);
 409
 410   SKIP: {
 411     $result = eval {
 412       is($db->delete_groups(1,2,3,4,5),5);
 413       my @features = $db->get_feature_by_name(Sequence => 'Contig2');
 414       is($db->delete_groups(@features),1);
 415       1;
 416     };
 417     if (!$result && $@ =~ /not implemented/i) {
 418       skip("delete_groups() not implemented by this adaptor",2);
 419     }
 420   }
 421
 422   # test ability to pass adaptors across a fork
 423   if (my $child = open(F,"-|")) { # parent reads from child
 424       ok(scalar <F>);
 425       close F;
 426   }
 427   else { # in child
 428       $db->clone;
 429       my @f = $db->features();
 430       print @f>0;
 431       exit 0;
 432   }
 433
 434   ok(!defined eval{$db->delete()});
 435   ok($db->delete(-force=>1));
 436   is(scalar $db->features,0);
 437   ok(!$db->segment('Contig1'));
 438
 439 }
 440
 441 }
 442
 443
 444
 445 END {
 446   unlink $fasta_files."/directory.index";
 447 }