t/LocalDB/BioDBGFF.t

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