t/BioDBSeqFeature.t

   1 #-*-Perl-*-
   2 ## Bioperl Test Harness Script for Modules
   3
   4 # Before `make install' is performed this script should be runnable with
   5 # `make test'. After `make install' it should work as `perl test.t'
   6
   7 use strict;
   8 use Bio::Root::IO;
   9 use FindBin '$Bin';
  10 use constant TEST_COUNT => 52;
  11 use constant GFF_FILE    => Bio::Root::IO->catfile('t','data',
  12                                            'seqfeaturedb','test.gff3');
  13
  14 BEGIN {
  15     # to handle systems with no installed Test module
  16     # we include the t dir (where a copy of Test.pm is located)
  17     # as a fallback
  18     eval { require Test; };
  19     if( $@ ) {
  20         use lib 't';
  21     }
  22     use Test;
  23     plan test => TEST_COUNT;
  24     $ENV{ORACLE_HOME} ||= '/home/oracle/Home';
  25 }
  26
  27 use lib '.','..','./blib/lib';
  28 use Bio::DB::SeqFeature::Store;
  29 use Bio::DB::SeqFeature::Store::GFF3Loader;
  30
  31 sub bail ($;$) {
  32   my $count = shift;
  33   my $explanation = shift;
  34   for (1..$count) {
  35     skip($explanation,1);
  36   }
  37   exit 0;
  38 }
  39
  40 sub fail ($) {
  41   my $count = shift;
  42   for (1..$count) {
  43     ok(0);
  44   }
  45   exit 0;
  46 }
  47
  48 my (@f,$f,@s,$s,$seq1,$seq2);
  49
  50 my @args = @ARGV;
  51 @args = (-adaptor => 'memory') unless @args;
  52
  53 my $db = eval { Bio::DB::SeqFeature::Store->new(@args) };
  54 warn $@ if $@;
  55 ok($db);
  56 fail(TEST_COUNT - 1) unless $db;
  57
  58 my $loader = eval { Bio::DB::SeqFeature::Store::GFF3Loader->new(-store=>$db) };
  59 warn $@ if $@;
  60 ok($loader);
  61 fail(TEST_COUNT - 2) unless $loader;
  62
  63 # exercise the loader
  64 ok($loader->load(GFF_FILE));
  65
  66 # there should be one gene named 'abc-1'
  67 @f = $db->get_features_by_name('abc-1');
  68 ok(@f==1);
  69
  70 $f = $f[0];
  71 # there should be three subfeatures of type "exon" and three of type "CDS"
  72 ok($f->get_SeqFeatures('exon')==3);
  73 ok($f->get_SeqFeatures('CDS')==3);
  74
  75 # the sequence of feature abc-1 should match the sequence of the first exon at the beginning
  76 $seq1 = $f->seq->seq;
  77 $seq2 = (sort {$a->start<=>$b->start} $f->get_SeqFeatures('exon'))[0]->seq->seq;
  78 ok(substr($seq1,0,length $seq2) eq $seq2);
  79
  80 # sequence lengths should match
  81 ok(length $seq1 == $f->length);
  82
  83 # if we pull out abc-1 again we should get the same object
  84 ($s) = $db->get_features_by_name('abc-1');
  85 ok($f eq $s);
  86
  87 # we should get two objects when we ask for abc-1 using get_features_by_alias
  88 # this also depends on selective subfeature indexing
  89 @f = $db->get_features_by_alias('abc-1');
  90 ok(@f==2);
  91
  92 # the two features should be different
  93 ok($f[0] ne $f[1]);
  94
  95 # test that targets are working
  96 ($f) = $db->get_features_by_name('match1');
  97 ok(defined $f);
  98 $s = $f->target;
  99 ok(defined $s);
 100 ok($s->seq_id  eq 'CEESC13F');
 101 $seq1 = $s->seq->seq;
 102 ok(substr($seq1,0,10) eq 'ttgcgttcgg');
 103
 104 # can we fetch subfeatures?
 105 # gene3.a has the Index=1 attribute, so we should fetch it
 106 ($f) = $db->get_features_by_name('gene3.a');
 107 ok($f);
 108
 109 # gene 3.b doesn't have an index, so we shouldn't get it
 110 ($f) = $db->get_features_by_name('gene3.b');
 111 ok(!$f);
 112
 113 # test three-tiered genes
 114 ($f) = $db->get_features_by_name('gene3');
 115 ok($f);
 116 my @transcripts = $f->get_SeqFeatures;
 117 ok(@transcripts == 2);
 118 ok($transcripts[0]->method eq 'mRNA');
 119 ok($transcripts[0]->source eq 'confirmed');
 120
 121 # test that exon #2 is shared between the two transcripts
 122 my @exons1      = $transcripts[0]->get_SeqFeatures('CDS');
 123 ok(@exons1 == 3);
 124 my @exons2      = $transcripts[1]->get_SeqFeatures('CDS');
 125 my ($shared1)   = grep {$_->display_name||'' eq 'shared_exon'} @exons1;
 126 my ($shared2)   = grep {$_->display_name||'' eq 'shared_exon'} @exons2;
 127 ok($shared1 && $shared2);
 128 ok($shared1 eq $shared2);
 129 ok($shared1->primary_id eq $shared2->primary_id);
 130
 131 # test attributes
 132 ok($shared1->phase == 0);
 133 ok($shared1->strand eq +1);
 134 ok(($f->attributes('expressed'))[0] eq 'yes');
 135
 136 # test autoloading
 137 my ($gene3a) = grep { $_->display_name eq 'gene3.a'} @transcripts;
 138 my ($gene3b) = grep { $_->display_name eq 'gene3.b'} @transcripts;
 139 ok($gene3a);
 140 ok($gene3b);
 141 ok($gene3a->Is_expressed);
 142 ok(!$gene3b->Is_expressed);
 143
 144 # the representation of the 3'-UTR in the two transcripts a and b is
 145 # different (not recommended but supported by the GFF3 spec). In the
 146 # first case, there are two 3'UTRs existing as independent
 147 # features. In the second, there is one UTR with a split location.
 148 ok($gene3a->Three_prime_UTR == 2);
 149 ok($gene3b->Three_prime_UTR == 1);
 150 my ($utr) = $gene3b->Three_prime_UTR;
 151 ok($utr->segments == 2);
 152 my $location = $utr->location;
 153 ok($location->isa('Bio::Location::Split'));
 154 ok($location->sub_Location == 2);
 155
 156 # ok, test that queries are working properly.
 157 # find all features with the attribute "expressed"
 158 @f = $db->get_features_by_attribute({expressed=>'yes'});
 159 ok(@f == 2);
 160
 161 # find all top-level features on Contig3 -- there should be two
 162 @f = $db->get_features_by_location(-seq_id=>'Contig3');
 163 ok(@f == 2);
 164
 165 # find all top-level features on Contig3 of type 'assembly_component'
 166 @f = $db->features(-seq_id=>'Contig3',-type=>'assembly_component');
 167 ok(@f==1);
 168
 169 # test iteration
 170 @f = $db->features;
 171 my $feature_count = @f;
 172 ok($feature_count > 0);
 173
 174 my $i = $db->get_seq_stream;
 175 ok($i);
 176
 177 my $count;
 178 while ($i->next_seq) { $count++ }
 179 ok($feature_count == $count);
 180
 181 # regression test on bug in which get_SeqFeatures('type') did not filter inline segments
 182 @f = $db->get_features_by_name('agt830.3');
 183 ok(@f && !$f[0]->get_SeqFeatures('exon'));
 184 ok(@f && $f[0]->get_SeqFeatures('EST_match'));
 185
 186 # regression test on bug in which the load_id disappeared
 187 ok(@f && $f[0]->load_id eq 'Match2');
 188
 189 # regress on proper handling of multiple ID features
 190 my ($alignment) = $db->get_features_by_name('agt830.5');
 191 ok($alignment);
 192 ok($alignment->target->start == 1 && $alignment->target->end == 654);
 193 ok($alignment->get_SeqFeatures == 2);
 194 my $gff3 = $alignment->gff3_string(1);
 195 my @lines = split "\n",$gff3;
 196 ok (@lines == 2);
 197 ok ("@lines" !~ /Parent=/s);
 198 ok ("@lines" =~ /ID=/s);
 199
 200 1;
 201
 202 __END__
 203