t/SeqFeature/Unflattener.t

   1 # -*-Perl-*- Test Harness script for Bioperl
   2 # $Id$
   3
   4 use strict;
   5
   6 BEGIN {
   7     use Bio::Root::Test;
   8
   9     test_begin(-tests => 21);
  10
  11     use_ok('Bio::SeqIO');
  12     use_ok('Bio::SeqFeature::Tools::Unflattener');
  13 }
  14
  15 my $verbosity = test_debug();
  16
  17 my ($seq, @sfs);
  18 my $unflattener = Bio::SeqFeature::Tools::Unflattener->new;
  19 $unflattener->verbose($verbosity);
  20
  21
  22 if (1) {
  23     my @path = ("NC_000007-ribosomal-slippage.gb");
  24     $seq = getseq(@path);
  25
  26     my @topsfs = $seq->get_SeqFeatures;
  27     if( $verbosity > 0 ) {
  28         warn sprintf "TOP:%d\n", scalar(@topsfs);
  29         write_hier(@topsfs);
  30     }
  31
  32     # UNFLATTEN
  33     @sfs = $unflattener->unflatten_seq(-seq       => $seq,
  34                                        -use_magic => 1);
  35     if( $verbosity > 0 ) {
  36         warn "\n\nPOST PROCESSING:\n";
  37         write_hier(@sfs);
  38         warn sprintf "PROCESSED:%d\n", scalar(@sfs);
  39     }
  40     is(@sfs, 1995);
  41 }
  42
  43
  44 if (1) {
  45     my @path = ("ribosome-slippage.gb");
  46     $seq     = getseq(@path);
  47
  48     my @topsfs = $seq->get_SeqFeatures;
  49     if( $verbosity > 0 ) {
  50         warn sprintf "TOP:%d\n", scalar(@topsfs);
  51         write_hier(@topsfs);
  52     }
  53
  54     # UNFLATTEN
  55     @sfs = $unflattener->unflatten_seq(-seq       => $seq,
  56                                        -use_magic => 1);
  57     if( $verbosity > 0 ) {
  58         warn "\n\nPOST PROCESSING:\n";
  59         write_hier(@sfs);
  60         warn sprintf "PROCESSED:%d\n", scalar(@sfs);
  61     }
  62     is(@sfs, 3);
  63 }
  64
  65
  66 if (1) {
  67     my @path = ("AE003644_Adh-genomic.gb");
  68     $seq     = getseq(@path);
  69
  70     is ($seq->accession_number, 'AE003644');
  71     my @topsfs = $seq->get_SeqFeatures;
  72     if( $verbosity > 0 ) {
  73         warn sprintf "TOP:%d\n", scalar(@topsfs);
  74         write_hier(@topsfs);
  75     }
  76
  77     # UNFLATTEN
  78     @sfs = $unflattener->unflatten_seq(-seq       => $seq,
  79                                        -group_tag => 'locus_tag');
  80     if( $verbosity > 0 ) {
  81         warn "\n\nPOST PROCESSING:\n";
  82         write_hier(@sfs);
  83         warn sprintf "PROCESSED:%d\n", scalar(@sfs);
  84     }
  85     is(@sfs, 21);
  86 }
  87
  88 # now try again, using a custom subroutine to link together features
  89 $seq = getseq("AE003644_Adh-genomic.gb");
  90 @sfs = $unflattener->unflatten_seq
  91     (-seq       => $seq,
  92      -group_tag => 'locus_tag',
  93      -resolver_method =>
  94         sub {
  95              my $self = shift;
  96              my ($sf, @candidate_container_sfs) = @_;
  97
  98              if ($sf->has_tag('note')) {
  99                  my @notes   = $sf->get_tag_values('note');
 100                  my @trnames = map {/from transcript\s+(.*)/; $1;}
 101                                @notes;
 102                  @trnames    = grep {$_} @trnames;
 103
 104                  my $trname;
 105                  if (@trnames == 0) {
 106                      $self->throw("UNRESOLVABLE");
 107                  }
 108                  elsif (@trnames == 1) {
 109                      $trname = $trnames[0];
 110                  }
 111                  else {
 112                      $self->throw("AMBIGUOUS: @trnames");
 113                  }
 114
 115                  my @container_sfs =
 116                      grep {
 117                            my ($product) =
 118                                  $_->has_tag('product') ? $_->get_tag_values('product')
 119                                : ('');
 120                            $product eq $trname;
 121                      } @candidate_container_sfs;
 122
 123                  if (@container_sfs == 0) {
 124                      $self->throw("UNRESOLVABLE");
 125                  }
 126                  elsif (@container_sfs == 1) {
 127                      # we got it!
 128                      return ($container_sfs[0]=>0);
 129                  }
 130                  else {
 131                      $self->throw("AMBIGUOUS");
 132                  }
 133              }
 134         });
 135 $unflattener->feature_from_splitloc(-seq => $seq);
 136 if( $verbosity > 0 ) {
 137     warn "\n\nPOST PROCESSING:\n";
 138     write_hier(@sfs);
 139     warn sprintf "PROCESSED2:%d\n", scalar(@sfs);
 140 }
 141 is(@sfs, 21);
 142
 143 # try again; different sequence
 144 # this is an E-Coli seq with no mRNA features;
 145 # we just want to link all features directly with gene
 146
 147 $seq = getseq("D10483.gbk");
 148
 149 # UNFLATTEN
 150 @sfs = $unflattener->unflatten_seq(-seq       => $seq,
 151                                    -partonomy => {'*'=>'gene'});
 152 if( $verbosity > 0 ) {
 153     warn "\n\nPOST PROCESSING:\n";
 154     write_hier(@sfs);
 155     warn sprintf "PROCESSED:%d\n", scalar(@sfs);
 156 }
 157 is(@sfs, 98);
 158
 159 # this sequence has no locus_tag or or gene tags
 160 $seq = getseq("AY763288.gb");
 161
 162 # UNFLATTEN
 163 @sfs = $unflattener->unflatten_seq(-seq       => $seq,
 164                                    -use_magic => 1);
 165 if( $verbosity > 0 ) {
 166     warn "\n\nPOST PROCESSING:\n";
 167     write_hier(@sfs);
 168     warn sprintf "PROCESSED:%d\n", scalar(@sfs);
 169 }
 170 is(@sfs, 3);
 171
 172 # try again; different sequence - dicistronic gene, mRNA record
 173
 174 $seq = getseq("X98338_Adh-mRNA.gb");
 175
 176 # UNFLATTEN
 177 @sfs = $unflattener->unflatten_seq(-seq       => $seq,
 178                                    -partonomy => {'*'=>'gene'});
 179 if( $verbosity > 0 ) {
 180     warn "\n\nPOST PROCESSING:\n";
 181     write_hier(@sfs);
 182     warn sprintf "PROCESSED:%d\n", scalar(@sfs);
 183 }
 184 is(@sfs, 7);
 185
 186 # try again; this sequence has no CDSs but rRNA present
 187
 188 $seq = getseq("no_cds_example.gb");
 189
 190 # UNFLATTEN
 191 @sfs = $unflattener->unflatten_seq(-seq       => $seq,
 192                                    -use_magic => 1);
 193 if( $verbosity > 0 ) {
 194     warn "\n\nPOST PROCESSING:\n";
 195     write_hier(@sfs);
 196     warn sprintf "PROCESSED:%d\n", scalar(@sfs);
 197 }
 198
 199 my @all_sfs = $seq->get_all_SeqFeatures;
 200
 201 my @exons = grep { $_-> primary_tag eq 'exon' } @all_sfs ;
 202
 203 is(@exons, 2);
 204
 205
 206 if (1) {
 207     # this is an arabidopsise gbk record. it has no mRNA features.
 208     # it has explicit exon/intron records
 209     my @path = ("ATF14F8.gbk");
 210     $seq     = getseq(@path);
 211
 212     is ($seq->accession_number, 'AL391144');
 213     my @topsfs = $seq->get_SeqFeatures;
 214     my @cdss   = grep {$_->primary_tag eq 'CDS'} @topsfs;
 215     my $n      = scalar(@topsfs);
 216     if( $verbosity > 0 ) {
 217         warn sprintf "TOP:%d\n", scalar(@topsfs);
 218         write_hier(@topsfs);
 219     }
 220
 221     # UNFLATTEN
 222     @sfs = $unflattener->unflatten_seq(-seq       => $seq,
 223                                        -use_magic => 1);
 224     @sfs = $seq->get_SeqFeatures;
 225     if( $verbosity > 0 ) {
 226         warn "\n\nPOST PROCESSING:\n";
 227         write_hier(@sfs);
 228         warn sprintf "PROCESSED/TOP:%d\n", scalar(@sfs);
 229     }
 230     is(@sfs,28);
 231
 232     my @allsfs = $seq->get_all_SeqFeatures;
 233     is(@allsfs,202);
 234
 235     my @mrnas = grep {$_->primary_tag eq 'mRNA'} @allsfs;
 236     if( $verbosity > 0 ) {
 237         warn sprintf "ALL:%d\n",   scalar(@allsfs);
 238         warn sprintf "mRNAs:%d\n", scalar(@mrnas);
 239     }
 240     # relationship between mRNA and CDS should be one-one
 241     is(@mrnas,@cdss);
 242 }
 243
 244
 245 if (1) {
 246     # this is a record from FlyBase
 247     # it has mRNA features, and explicit exon/intron records
 248     my @path = ("AnnIX-v003.gbk");
 249     $seq     = getseq(@path);
 250
 251     my @topsfs = $seq->get_SeqFeatures;
 252     if( $verbosity > 0 ) {
 253         warn sprintf "TOP:%d\n", scalar(@topsfs);
 254         write_hier(@topsfs);
 255     }
 256
 257     # UNFLATTEN
 258     @sfs = $unflattener->unflatten_seq(-seq       => $seq,
 259                                        -use_magic => 1);
 260     @sfs = $seq->get_SeqFeatures;
 261     if( $verbosity > 0 ) {
 262         warn "\n\nPOST PROCESSING:\n";
 263         write_hier(@sfs);
 264         warn sprintf "PROCESSED/TOP:%d\n", scalar(@sfs);
 265     }
 266     is scalar(@sfs), 1;
 267
 268     my @exons = grep {$_->primary_tag eq 'exon'} $seq->get_all_SeqFeatures;
 269     is scalar(@exons), 6;    # total number of exons per splice
 270
 271     my %numberh = map {$_->get_tag_values("number") => 1} @exons;
 272     my @numbers = keys %numberh;
 273     if( $verbosity > 0 ) {
 274         warn sprintf "DISTINCT EXONS: %d [@numbers]\n", scalar(@numbers);
 275     }
 276     is scalar(@numbers), 6;  # distinct exons
 277 }
 278
 279
 280 if (1) {
 281     # example of a BAD genbank entry
 282     my @path = ("dmel_2Lchunk.gb");
 283     $seq     = getseq(@path);
 284
 285     my @topsfs = $seq->get_SeqFeatures;
 286     if( $verbosity > 0 ) {
 287         warn sprintf "TOP:%d\n", scalar(@topsfs);
 288         write_hier(@topsfs);
 289     }
 290
 291     # UNFLATTEN
 292     #
 293     # we EXPECT problems with this erroneous record
 294     $unflattener->error_threshold(2);
 295     @sfs = $unflattener->unflatten_seq(-seq       => $seq,
 296                                        -use_magic => 1);
 297     @sfs = $seq->get_SeqFeatures;
 298     if( $verbosity > 0 ) {
 299         warn "\n\nPOST PROCESSING:\n";
 300         write_hier(@sfs);
 301         warn sprintf "PROCESSED/TOP:%d\n", scalar(@sfs);
 302     }
 303     is scalar(@sfs), 2;
 304
 305     my @exons = grep {$_->primary_tag eq 'exon'} $seq->get_all_SeqFeatures;
 306     is scalar(@exons), 2;    # total number of exons per splice
 307
 308     my @probs = $unflattener->get_problems;
 309     $unflattener->report_problems(\*STDERR) if $verbosity > 0;
 310     $unflattener->clear_problems;
 311     if( $verbosity > 0 ) {
 312         warn sprintf "PROBLEMS ENCOUNTERED: %d (EXPECTED: 6)\n", scalar(@probs);
 313     }
 314     is scalar(@probs), 6;
 315 }
 316
 317
 318 sub write_hier {
 319     my @sfs = @_;
 320     _write_hier(0, @sfs);
 321 }
 322
 323 sub _write_hier {
 324     my $indent = shift;
 325     my @sfs = @_;
 326     foreach my $sf (@sfs) {
 327         my $label = '?';
 328         if ($sf->has_tag('gene')) {
 329             ($label) = $sf->get_tag_values('gene');
 330         }
 331         elsif ($sf->has_tag('product')) {
 332             ($label) = $sf->get_tag_values('product');
 333         }
 334         elsif ($sf->has_tag('number')) {
 335             $label = join("; ", $sf->get_tag_values('number'));
 336         }
 337         warn sprintf "%s%s $label\n", '  ' x $indent, $sf->primary_tag;
 338         my @sub_sfs = $sf->sub_SeqFeature;
 339         _write_hier($indent+1, @sub_sfs);
 340     }
 341 }
 342
 343 sub getseq {
 344     my @path  = @_;
 345     my $seqio = Bio::SeqIO->new('-file'   => test_input_file(@path),
 346                                 '-format' => 'GenBank');
 347     $seqio->verbose($verbosity);
 348
 349     my $seq = $seqio->next_seq();
 350     return $seq;
 351 }