t/Unflattener2.t

   1 # -*-Perl-*- mode (to keep my emacs happy)
   2 # $Id$
   3
   4 use strict;
   5 use vars qw($DEBUG $TESTCOUNT);
   6 BEGIN {
   7     eval { require Test; };
   8     if( $@ ) {
   9         use lib 't';
  10     }
  11     use Test;
  12     $TESTCOUNT = 11;
  13     plan tests => $TESTCOUNT;
  14 }
  15
  16 use Bio::Seq;
  17 use Bio::SeqIO;
  18 use Bio::Root::IO;
  19 use Bio::SeqFeature::Tools::Unflattener;
  20
  21 ok(1);
  22
  23 my $verbosity = -1;   # Set to -1 for release version, so warnings aren't printed
  24 #$verbosity = 1;
  25
  26 my ($seq, @sfs);
  27 my $unflattener = Bio::SeqFeature::Tools::Unflattener->new;
  28 $unflattener->verbose($verbosity);
  29
  30 if (1) {
  31
  32     # this is an arabidopsise gbk record. it has no mRNA features.
  33     # it has explicit exon/intron records
  34
  35     my @path = ("t","data","ATF14F8.gbk");
  36     $seq = getseq(@path);
  37
  38     ok ($seq->accession_number, 'AL391144');
  39     my @topsfs = $seq->get_SeqFeatures;
  40     my @cdss = grep {$_->primary_tag eq 'CDS'} @topsfs;
  41     my $n = scalar(@topsfs);
  42     if( $verbosity > 0 ) {
  43         warn sprintf "TOP:%d\n", scalar(@topsfs);
  44         write_hier(@topsfs);
  45     }
  46     # UNFLATTEN
  47     @sfs = $unflattener->unflatten_seq(-seq=>$seq,
  48                                        -use_magic=>1,
  49                                       );
  50     @sfs = $seq->get_SeqFeatures;
  51     if( $verbosity > 0 ) {
  52         warn "\n\nPOST PROCESSING:\n";
  53         write_hier(@sfs);
  54         warn sprintf "PROCESSED/TOP:%d\n", scalar(@sfs);
  55     }
  56     ok(@sfs == 28);
  57     my @allsfs = $seq->get_all_SeqFeatures;
  58     ok(@allsfs == 202);
  59     my @mrnas = grep {$_->primary_tag eq 'mRNA'} @allsfs;
  60     if( $verbosity > 0 ) {
  61         warn sprintf "ALL:%d\n", scalar(@allsfs);
  62         warn sprintf "mRNAs:%d\n", scalar(@mrnas);
  63     }
  64
  65     # relationship between mRNA and CDS should be one-one
  66     ok(@mrnas == @cdss);
  67 }
  68
  69 if (1) {
  70
  71     # this is a record from FlyBase
  72     # it has mRNA features, and explicit exon/intron records
  73
  74     my @path = ("t","data","AnnIX-v003.gbk");
  75     $seq = getseq(@path);
  76
  77     my @topsfs = $seq->get_SeqFeatures;
  78     if( $verbosity > 0 ) {
  79         warn sprintf "TOP:%d\n", scalar(@topsfs);
  80         write_hier(@topsfs);
  81     }
  82     # UNFLATTEN
  83     @sfs = $unflattener->unflatten_seq(-seq=>$seq,
  84                                        -use_magic=>1,
  85                                       );
  86     @sfs = $seq->get_SeqFeatures;
  87     if( $verbosity > 0 ) {
  88         warn "\n\nPOST PROCESSING:\n";
  89         write_hier(@sfs);
  90         warn sprintf "PROCESSED/TOP:%d\n", scalar(@sfs);
  91     }
  92     ok scalar(@sfs), 1;
  93     my @exons = grep {$_->primary_tag eq 'exon'} $seq->get_all_SeqFeatures;
  94     ok scalar(@exons), 6;    # total number of exons per splice
  95     my %numberh = map {$_->get_tag_values("number") => 1} @exons;
  96     my @numbers = keys %numberh;
  97     if( $verbosity > 0 ) {
  98         warn sprintf "DISTINCT EXONS: %d [@numbers]\n", scalar(@numbers);
  99     }
 100     ok scalar(@numbers), 6;  # distinct exons
 101 }
 102
 103 if (1) {
 104
 105     # example of a BAD genbank entry
 106
 107     my @path = ("t","data","dmel_2Lchunk.gb");
 108     $seq = getseq(@path);
 109
 110     my @topsfs = $seq->get_SeqFeatures;
 111     if( $verbosity > 0 ) {
 112         warn sprintf "TOP:%d\n", scalar(@topsfs);
 113         write_hier(@topsfs);
 114     }
 115     # UNFLATTEN
 116     #
 117     # we EXPECT problems with this erroneous record
 118     $unflattener->error_threshold(2);
 119     @sfs = $unflattener->unflatten_seq(-seq=>$seq,
 120                                        -use_magic=>1,
 121                                       );
 122     my @probs = $unflattener->get_problems;
 123     $unflattener->report_problems(\*STDERR) if $verbosity > 0;
 124     $unflattener->clear_problems;
 125     @sfs = $seq->get_SeqFeatures;
 126     if( $verbosity > 0 ) {
 127         warn "\n\nPOST PROCESSING:\n";
 128         write_hier(@sfs);
 129         warn sprintf "PROCESSED/TOP:%d\n", scalar(@sfs);
 130     }
 131     ok scalar(@sfs), 2;
 132     my @exons = grep {$_->primary_tag eq 'exon'} $seq->get_all_SeqFeatures;
 133     ok scalar(@exons), 2;    # total number of exons per splice
 134     if( $verbosity > 0 ) {
 135         warn sprintf "PROBLEMS ENCOUNTERED: %d (EXPECTED: 6)\n", scalar(@probs);
 136     }
 137     ok scalar(@probs), 6;
 138 }
 139
 140
 141 sub write_hier {
 142     my @sfs = @_;
 143     _write_hier(0, @sfs);
 144 }
 145
 146 sub _write_hier {
 147     my $indent = shift;
 148     my @sfs = @_;
 149     foreach my $sf (@sfs) {
 150         my $label = '?';
 151         if ($sf->has_tag('gene')) {
 152             ($label) = $sf->get_tag_values('gene');
 153         }
 154         if ($sf->has_tag('product')) {
 155             ($label) = $sf->get_tag_values('product');
 156         }
 157         if ($sf->has_tag('number')) {
 158             $label = join("; ", $sf->get_tag_values('number'));
 159         }
 160         printf "%s%s $label\n", '  ' x $indent, $sf->primary_tag;
 161         my @sub_sfs = $sf->sub_SeqFeature;
 162         _write_hier($indent+1, @sub_sfs);
 163     }
 164 }
 165
 166 sub getseq {
 167     my @path = @_;
 168     my $seqio =
 169       Bio::SeqIO->new('-file'=> Bio::Root::IO->catfile(
 170                                                        @path
 171                                                       ),
 172                       '-format' => 'GenBank');
 173     $seqio->verbose($verbosity);
 174
 175     my $seq = $seqio->next_seq();
 176     return $seq;
 177 }
 178
 179 # 1 2,3
 180 # 2 1,2
 181 # 3 4,5
 182 # 4 1,4,5,6
 183 # 5 1,4,5,6
 184 # 6 1,4,5,6