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