t/BioDBSeqFeature.t

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