t/SeqFeature/Unflattener.t

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