← Back to branch summary

~ubuntu-branches/ubuntu/raring/bioperl/raring

« back to all changes in this revision

Viewing changes to Bio/SeqFeature/Tools/Unflattener.pm

  • Committer: Bazaar Package Importer
  • Author(s): Charles Plessy
  • Date: 2008-03-18 14:44:57 UTC
  • mfrom: (4 hardy)
  • mto: This revision was merged to the branch mainline in revision 6.
  • Revision ID: james.westby@ubuntu.com-20080318144457-1jjoztrvqwf0gruk
http://bugs.debian.org/448890
* debian/control:
  - Removed MIA Matt Hope (dopey) from the Uploaders field.
    Thank you for your work, Matt. I hope you are doing well.
  - Downgraded some recommended package to the 'Suggests' priority,
    according to the following discussion on Upstream's mail list.
    http://bioperl.org/pipermail/bioperl-l/2008-March/027379.html
    (Closes: #448890)
* debian/copyright converted to machine-readable format.

Show diffs side-by-side

added added

removed removed

Lines of Context:
Bio/SeqFeature/Tools/Unflattener.pm
1
 
# $Id: Unflattener.pm,v 1.19 2003/12/16 22:31:16 cjm Exp $
 
1
# $Id: Unflattener.pm,v 1.37.4.2 2006/10/02 23:10:28 sendu Exp $
2
2
#
3
3
# bioperl module for Bio::SeqFeature::Tools::Unflattener
4
4
#
583
583
Bioperl modules. Send your comments and suggestions preferably to the
584
584
Bioperl mailing lists  Your participation is much appreciated.
585
585
 
586
 
  bioperl-l@bioperl.org                         - General discussion
587
 
  http://bio.perl.org/MailList.html             - About the mailing lists
 
586
  bioperl-l@bioperl.org                  - General discussion
 
587
  http://bioperl.org/wiki/Mailing_lists  - About the mailing lists
588
588
 
589
589
=head2 Reporting Bugs
590
590
 
591
591
report bugs to the Bioperl bug tracking system to help us keep track
592
 
 the bugs and their resolution.  Bug reports can be submitted via
593
 
 email or the web:
 
592
the bugs and their resolution.  Bug reports can be submitted via the
 
593
web:
594
594
 
595
 
  bioperl-bugs@bio.perl.org
596
 
  http://bugzilla.bioperl.org/
 
595
  http://bugzilla.open-bio.org/
597
596
 
598
597
=head1 AUTHOR - Chris Mungall
599
598
 
610
609
# Let the code begin...
611
610
 
612
611
package Bio::SeqFeature::Tools::Unflattener;
613
 
use vars qw(@ISA);
614
612
use strict;
615
613
 
616
614
# Object preamble - inherits from Bio::Root::Root
617
 
use Bio::Root::Root;
618
615
use Bio::Location::Simple;
619
616
use Bio::SeqFeature::Generic;
620
617
use Bio::Range;
621
618
 
622
619
 
623
 
@ISA = qw(Bio::Root::Root);
 
620
use base qw(Bio::Root::Root);
624
621
 
625
622
=head2 new
626
623
 
666
663
    my @probs = $self->get_problems;
667
664
    if (!$self->{_problems_reported} &&
668
665
        scalar(@probs)) {
669
 
        print STDERR 
670
 
          "WARNING: There are UNREPORTED PROBLEMS.\n".
 
666
        $self->warn(
 
667
            "WARNING: There are UNREPORTED PROBLEMS.\n".
671
668
            "You may wish to use the method report_problems(), \n",
672
 
              "or ignore_problems() on the Unflattener object\n";
 
669
            "or ignore_problems() on the Unflattener object\n");
673
670
    }
674
671
    return;
675
672
}
774
771
            CDS => 'mRNA',
775
772
            exon => 'mRNA',
776
773
            intron => 'mRNA',
 
774
 
 
775
            pseudoexon => 'pseudogene',
 
776
            pseudointron => 'pseudogene',
 
777
            pseudotranscript => 'pseudogene',
777
778
           };
778
779
}
779
780
 
828
829
  source
829
830
  gene
830
831
  CDS
831
 
  exon
832
 
  intron
 
832
  exon [optional]
 
833
  intron [optional]
833
834
 
834
835
there are no mRNA features
835
836
 
841
842
      intron
842
843
 
843
844
exon and intron may or may not be present; they may be implicit from
844
 
the mRNA 'join' location
 
845
the CDS 'join' location
845
846
 
846
847
=back
847
848
 
907
908
    my $self = shift;
908
909
 
909
910
    $self->{'_problems'} = [] unless exists($self->{'_problems'});
910
 
    if ($self->verbose) {
911
 
        print "PROBLEM: $_\n" foreach @_;
 
911
    if ($self->verbose > 0) {
 
912
        warn( "PROBLEM: $_\n") foreach @_;
912
913
    }
913
914
    push(@{$self->{'_problems'}}, @_);
914
915
}
1143
1144
       # use magic to guess the group tag
1144
1145
       my @sfs_with_locus_tag =
1145
1146
         grep {$_->has_tag("locus_tag")} @flat_seq_features;
 
1147
       my @sfs_with_gene_tag =
 
1148
         grep {$_->has_tag("gene")} @flat_seq_features;
 
1149
       my @sfs_with_product_tag =
 
1150
         grep {$_->has_tag("product")} @flat_seq_features;
1146
1151
       if (@sfs_with_locus_tag) {
1147
1152
           if ($group_tag && $group_tag ne 'locus_tag') {
1148
1153
               $self->throw("You have explicitly set group_tag to be '$group_tag'\n".
1151
1156
                            "You can resolve this by either NOT setting -group_tag\n".
1152
1157
                            "OR you can unset -use_magic to regain control");
1153
1158
           }
1154
 
           # I know what's best for you.
1155
 
           # You want to be using /locus_tag=foo here
1156
 
           #
 
1159
 
 
1160
           # use /locus_tag instead of /gene tag for grouping
1157
1161
           # see GenBank entry AE003677 (version 3) for an example
1158
1162
           $group_tag = 'locus_tag';
1159
 
           if ($self->verbose) {
1160
 
               print "Set group tag to: $group_tag\n";
1161
 
           }
 
1163
           if ($self->verbose > 0) {
 
1164
               warn "Set group tag to: $group_tag\n";
 
1165
           }
 
1166
       }
 
1167
 
 
1168
       # on rare occasions, records will have no /gene or /locus_tag
 
1169
       # but it WILL have /product tags. These serve the same purpose
 
1170
       # for grouping. For an example, see AY763288 (also in t/data)
 
1171
       if (@sfs_with_locus_tag==0 &&
 
1172
           @sfs_with_gene_tag==0 &&
 
1173
           @sfs_with_product_tag>0 &&
 
1174
           !$group_tag) {
 
1175
           $group_tag = 'product';
 
1176
           if ($self->verbose > 0) {
 
1177
               warn "Set group tag to: $group_tag\n";
 
1178
           }
 
1179
           
1162
1180
       }
1163
1181
   }
1164
1182
   if (!$group_tag) {
1257
1275
 
1258
1276
   # -
1259
1277
 
 
1278
   # PSEUDOGENES, PSEUDOEXONS AND PSEUDOINTRONS
 
1279
   # these are indicated with the /pseudo tag
 
1280
   # these are mapped to a different type; they should NOT
 
1281
   # be treated as normal genes
 
1282
   foreach my $sf (@all_seq_features) {
 
1283
       if ($sf->has_tag('pseudo')) {
 
1284
           my $type = $sf->primary_tag;
 
1285
           # SO type is typically the same as the normal
 
1286
           # type but preceeded by "pseudo"
 
1287
           if ($type eq 'misc_RNA') {
 
1288
               $sf->primary_tag("pseudotranscript");
 
1289
           }
 
1290
           else {
 
1291
               $sf->primary_tag("pseudo$type");
 
1292
           }
 
1293
       }
 
1294
   }
 
1295
   # now some of the post-processing that follows which applies to
 
1296
   # genes will NOT be applied to pseudogenes; this is deliberate
 
1297
   # for example, gene models are normalised to be gene-transcript-exon
 
1298
   # for pseudogenes we leave them as pseudogene-pseudoexon
 
1299
 
1260
1300
   # --- MAGIC ---
1261
1301
   my $need_to_infer_exons = 0;
1262
1302
   my $need_to_infer_mRNAs = 0;
1276
1316
                        $_->has_tag($group_tag)} @flat_seq_features);
1277
1317
       my $n_cdss =
1278
1318
         scalar(grep {$_->primary_tag eq 'CDS'} @flat_seq_features);
1279
 
           
 
1319
       my $n_rnas =
 
1320
         scalar(grep {$_->primary_tag =~ /RNA/} @flat_seq_features);  
1280
1321
       # Are there any CDS features in the record?
1281
 
       if ($n_cdss) {
 
1322
       if ($n_cdss > 0) {
1282
1323
           # YES
1283
1324
           
1284
1325
           # - a pc gene model should contain at the least a CDS
1285
1326
 
1286
1327
           # Are there any mRNA features in the record?
1287
 
           if (!$n_mrnas) {
 
1328
           if ($n_mrnas == 0) {
 
1329
               # NO mRNAs:
1288
1330
               # looks like structure_type == 1
1289
1331
               $structure_type = 1;
1290
1332
               $need_to_infer_mRNAs = 1;
1291
1333
           }
1292
 
           elsif (!$n_mrnas_attached_to_gene) {
1293
 
               # The files _does_ have at least one mRNA feature,
1294
 
               # but none of them are part of a group, i.e. they
 
1334
           elsif ($n_mrnas_attached_to_gene == 0) {
 
1335
               # $n_mrnas > 0
 
1336
               # $n_mrnas_attached_to_gene = 0
 
1337
               #
 
1338
               # The entries _do_ contain mRNA features,
 
1339
               # but none of them are part of a group/gene, i.e. they
1295
1340
               # are 'floating'
1296
1341
 
1297
1342
               # this is an annoying weird file that has some floating
1298
1343
               # mRNA features; 
1299
1344
               # eg ftp.ncbi.nih.gov/genomes/Schizosaccharomyces_pombe/
 
1345
               
 
1346
               if ($self->verbose) {
 
1347
                   my @floating_mrnas =
 
1348
                     grep {$_->primary_tag eq 'mRNA' &&
 
1349
                             !$_->has_tag($group_tag)} @flat_seq_features;
 
1350
                   printf STDERR "Unattached mRNAs:\n";
 
1351
                   foreach my $mrna (@floating_mrnas) {
 
1352
                       $self->_write_sf_detail($mrna);
 
1353
                   }
 
1354
                   printf STDERR "Don't know how to deal with these; filter at source?\n";
 
1355
               }
1300
1356
 
1301
1357
               foreach (@flat_seq_features) {
1302
1358
                   if ($_->primary_tag eq 'mRNA') {
1319
1375
       }
1320
1376
       else {
1321
1377
           # this doesn't seem to be any kind of protein coding gene model
 
1378
           if ( $n_rnas > 0 ) {
 
1379
               $need_to_infer_exons = 1;
 
1380
           }
1322
1381
       }
1323
1382
 
1324
1383
       if ($need_to_infer_exons) {
1349
1408
   }
1350
1409
 
1351
1410
   # LOGGING
1352
 
   if ($self->verbose) {
1353
 
       print "GROUPS:\n";
 
1411
   if ($self->verbose > 0) {
 
1412
       printf STDERR "GROUPS:\n";
1354
1413
       foreach my $group (@groups) {
1355
1414
           $self->_write_group($group, $group_tag);
1356
1415
       }
1357
1416
   }
1358
1417
   # -
1359
1418
 
1360
 
   # >>>>>>>>>                   <<<<<<<<<<<<<
1361
1419
   # --------- FINISHED GROUPING -------------
1362
 
   # >>>>>>>>>                   <<<<<<<<<<<<<
1363
 
 
 
1420
 
 
1421
 
 
1422
   # TYPE CONTAINMENT HIERARCHY (aka partonomy)
1364
1423
   # set the containment hierarchy if desired
1365
1424
   # see docs for structure_type() method
1366
1425
   if ($structure_type) {
1367
1426
       if ($structure_type == 1) {
1368
1427
           $self->partonomy(
1369
 
                                        {CDS => 'gene',
1370
 
                                         exon => 'CDS',
1371
 
                                         intron => 'CDS',
1372
 
                                        }
1373
 
                                       );
 
1428
                            {CDS => 'gene',
 
1429
                             exon => 'CDS',
 
1430
                             intron => 'CDS',
 
1431
                            }
 
1432
                           );
1374
1433
       }
1375
1434
       else {
1376
1435
           $self->throw("structure_type $structure_type is currently unknown");
1387
1446
   }
1388
1447
 
1389
1448
   if ($use_magic) {
1390
 
       my @roots = $self->_get_partonomy_roots;
1391
1449
       # point all feature types without a container type to the root type.
1392
1450
       #
1393
1451
       # for example, if we have an unanticipated feature_type, say
1399
1457
                   my $type = $sf->primary_tag;
1400
1458
                   next if $type eq 'gene';
1401
1459
                   my $container_type = $self->get_container_type($type);
1402
 
                   my $root = $roots[0];
1403
1460
                   if (!$container_type) {
1404
 
                       $self->partonomy->{$type} = $root;
 
1461
                       $self->partonomy->{$type} = 'gene';
1405
1462
                   }
1406
1463
               }
1407
1464
           }
1435
1492
   $seq->remove_SeqFeatures;
1436
1493
   $seq->add_SeqFeature(@top_sfs);
1437
1494
 
1438
 
   # >>>>>>>>>                       <<<<<<<<<<<<<
1439
1495
   # --------- FINISHED UNFLATTENING -------------
1440
 
   # >>>>>>>>>                       <<<<<<<<<<<<<
1441
1496
 
1442
1497
   # lets see if there are any post-unflattening tasks we need to do
1443
1498
 
 
1499
   
 
1500
 
1444
1501
   # INFERRING mRNAs
1445
1502
   if ($need_to_infer_mRNAs) {
1446
 
       if ($self->verbose) {
1447
 
           print "** INFERRING mRNA from CDS\n";
 
1503
       if ($self->verbose > 0) {
 
1504
           printf STDERR "** INFERRING mRNA from CDS\n";
1448
1505
       }
1449
1506
       $self->infer_mRNA_from_CDS(-seq=>$seq);
1450
1507
   }
1517
1574
                           $exon->add_tag_value($tag, @vals);
1518
1575
                       }
1519
1576
                   }
1520
 
               } else {
 
1577
               } 
 
1578
               else {
 
1579
                   # no exons inferred at $locstr
1521
1580
                   push(@problems,
1522
1581
                        [1, 
1523
1582
                         "there is a conflict with exons; there was an explicitly ".
1541
1600
               my $thresh = $self->error_threshold;
1542
1601
               my @bad_problems = grep {$_->[0] > $thresh} @problems;
1543
1602
               if (@bad_problems) {
1544
 
                   print STDERR "PROBLEM:\n";
 
1603
                   printf STDERR "PROBLEM:\n";
1545
1604
                   $self->_write_hier(\@top_sfs);
1546
1605
                   # TODO - allow more fine grained control over this
1547
1606
                   $self->{_problems_reported} = 1;
1588
1647
    else {
1589
1648
        # @ranges > 1
1590
1649
        # split the group into disconnected ranges
1591
 
        if ($self->verbose) {
1592
 
            print "GROUP PRE-SPLIT:\n";
 
1650
        if ($self->verbose > 0) {
 
1651
            printf STDERR "GROUP PRE-SPLIT:\n";
1593
1652
            $self->_write_group($group, $self->group_tag);
1594
1653
        }
1595
1654
        @groups =
1599
1658
                  $_->intersection($range);
1600
1659
              } @sfs]
1601
1660
          } @ranges;
1602
 
        if ($self->verbose) {
1603
 
            print "SPLIT GROUPS:\n";
 
1661
        if ($self->verbose > 0) {
 
1662
            printf STDERR "SPLIT GROUPS:\n";
1604
1663
            $self->_write_group($_, $self->group_tag) foreach @groups;      
1605
1664
        }
1606
1665
    }
1639
1698
        # the latter is redundant with the CDS entry. So we shall get rid of
1640
1699
        # the latter with the following filter
1641
1700
 
1642
 
        if ($self->verbose) {
1643
 
            print "REMOVING DUPLICATES:\n";
 
1701
        if ($self->verbose > 0) {
 
1702
            printf STDERR "REMOVING DUPLICATES:\n";
1644
1703
        }
1645
1704
 
1646
1705
        @genes =
1756
1815
                          )],
1757
1816
                          @args);
1758
1817
 
1759
 
   if ($self->verbose) {
1760
 
       print "UNFLATTENING GROUP:\n";
 
1818
   if ($self->verbose > 0) {
 
1819
       printf STDERR "UNFLATTENING GROUP:\n";
1761
1820
       $self->_write_group($group, $self->group_tag);
1762
1821
   }
1763
1822
 
1801
1860
                 @container_sfs = 
1802
1861
                   grep {
1803
1862
                       my $match = 0;
1804
 
                       $self->_write_sf($_);
 
1863
                       $self->_write_sf($_) if $self->verbose > 0;
1805
1864
                       foreach my $tag (qw(product symbol label)) {
1806
1865
                           if ($_->has_tag($tag)) {
1807
1866
                               my @vals =
1843
1902
   # CONDITION: there must be at most one root
1844
1903
   if (@top_sfs > 1) {
1845
1904
       $self->_write_group($group, $self->group_tag);
1846
 
       print "TOP SFS:\n";
 
1905
       printf STDERR "TOP SFS:\n";
1847
1906
       $self->_write_sf($_) foreach @top_sfs;
1848
1907
       $self->throw("multiple top-sfs in group");
1849
1908
   }
1993
2052
       if (%unresolved) {
1994
2053
           my %childh = map {$_=>1} keys %unresolved;
1995
2054
           my %parenth = map {$_->[0]=>1} map {@$_} values %unresolved;
1996
 
           if ($self->verbose) {
1997
 
               printf "MATCHING %d CHILDREN TO %d PARENTS\n",
 
2055
           if ($self->verbose > 0) {
 
2056
               printf STDERR "MATCHING %d CHILDREN TO %d PARENTS\n",
1998
2057
                 scalar(keys %childh), scalar(keys %parenth);
1999
2058
           }
2000
2059
           # 99.99% of the time in genbank genomic record of structure type 0, we
2014
2073
       }
2015
2074
   }
2016
2075
 
2017
 
   # CONDITION:
2018
 
   # The graph %container 
2019
 
 
2020
2076
   # DEBUGGING CODE
2021
 
   if ($self->verbose && scalar(keys %unresolved)) {
2022
 
       print "UNRESOLVED PAIRS:\n";
 
2077
   if ($self->verbose > 0 && scalar(keys %unresolved)) {
 
2078
       printf STDERR "UNRESOLVED PAIRS:\n";
2023
2079
       foreach my $childsf (keys %unresolved) {
2024
2080
           my @poss = @{$unresolved{$childsf}};
2025
2081
           foreach my $p (@poss) {
2026
2082
               my $parentsf = $p->[0];
2027
2083
               $childsf = $idxsf{$childsf};
2028
 
               my @clabels = $childsf->get_tagset_values(qw(protein_id label product));
2029
 
               my @plabels = $parentsf->get_tagset_values(qw(transcript_id label product));
2030
 
               printf("  PAIR: $clabels[0] => $plabels[0]  (of %d)\n", 
2031
 
                      scalar(@poss));
 
2084
               my @clabels = ($childsf->get_tagset_values(qw(protein_id label product)), "?");
 
2085
               my @plabels = ($parentsf->get_tagset_values(qw(transcript_id label product)), "?");
 
2086
               printf STDERR
 
2087
                      ("  PAIR: $clabels[0] => $plabels[0]  (of %d)\n", 
 
2088
                       scalar(@poss));
2032
2089
           }
2033
2090
       }
2034
2091
   } # -- end of verbose
2053
2110
           $unresolved_problem_reported = 1;
2054
2111
       }
2055
2112
       foreach my $pair (@$new_pairs) {
2056
 
           if ($self->verbose) {
2057
 
               printf "  resolved pair @$pair\n";
 
2113
           if ($self->verbose > 0) {
 
2114
               printf STDERR "  resolved pair @$pair\n";
2058
2115
           }
2059
2116
           $container{$pair->[0]} = $pair->[1];
2060
2117
           delete $unresolved{$pair->[0]};
2073
2130
   foreach my $sf (@sfs) {
2074
2131
       my $container_sf = $container{$sf};
2075
2132
       if ($container_sf) {
2076
 
           eval {
2077
 
               # make $sf nested inside $container_sf
2078
 
               #
2079
 
               # this will throw an exception if $sf is
2080
 
               # not spatially within $container_sf
 
2133
           # make $sf nested inside $container_sf
 
2134
 
 
2135
           # first check if the container spatially contains the containee
 
2136
           if ($container_sf->contains($sf)) {
 
2137
               # add containee
2081
2138
               $container_sf->add_SeqFeature($sf);
2082
 
           };
2083
 
           if ($@) {
2084
 
               # probably a spatial containment problem
2085
 
               $self->problem(2,
2086
 
                              "bioperl add_SeqFeature says:$@",
 
2139
           }
 
2140
           else {
 
2141
               # weird case - the container does NOT spatially
 
2142
               # contain the containee;
 
2143
               # we expand and throw a warning
 
2144
               #
 
2145
               # for an example of this see ZFP91-CNTF dicistronic gene
 
2146
               # in NCBI chrom 11 build 34.3
 
2147
               $self->problem(1,
 
2148
                              "Container feature does not spatially contain ".
 
2149
                              "subfeature. Perhaps this is a dicistronic gene? ".
 
2150
                              "I am expanding the parent feature",
2087
2151
                              $container_sf,
2088
2152
                              $sf);
2089
 
           }
 
2153
               $container_sf->add_SeqFeature($sf, 'EXPAND');
 
2154
           }
2090
2155
       }
2091
2156
       else {
2092
2157
           push(@top, $sf);
2131
2196
 
2132
2197
    my $verbose = $self->verbose;
2133
2198
    #################################print "I";
2134
 
    if ($verbose) {
2135
 
        printf "find_best_matches: (/%d)\n", scalar(@$pairs);
 
2199
    if ($verbose > 0) {
 
2200
        printf STDERR "find_best_matches: (/%d)\n", scalar(@$pairs);
2136
2201
    }
2137
2202
 
2138
2203
    my %selected_children = map {($_->[0]=>1)} @$pairs;
2143
2208
    my %unresolved_parents = ();
2144
2209
    my %unresolved =
2145
2210
      map {
2146
 
          if ($verbose) {
2147
 
              printf "  $_ : %s\n", join("; ", map {"[@$_]"} @{$matrix->{$_}});
 
2211
          if ($verbose > 0) {
 
2212
              printf STDERR "  $_ : %s\n", join("; ", map {"[@$_]"} @{$matrix->{$_}});
2148
2213
          }
2149
2214
          if ($selected_children{$_}) {
2150
2215
              ();
2231
2296
    my $group_tag = shift || 'gene';
2232
2297
 
2233
2298
    my $f = $group->[0];
2234
 
    my $label = '';
 
2299
    my $label = '?';
2235
2300
    if ($f->has_tag($group_tag)) {
2236
2301
        ($label) = $f->get_tag_values($group_tag);
2237
2302
    }
2238
 
    printf("  GROUP [%s]:%s\n",
2239
 
           $label,
2240
 
           join(' ',
2241
 
                map { $_->primary_tag } @$group));
 
2303
    if( $self->verbose > 0 ) { 
 
2304
        printf STDERR ("  GROUP [%s]:%s\n",
 
2305
               $label,
 
2306
               join(' ',
 
2307
                    map { $_->primary_tag } @$group));
 
2308
    }
2242
2309
 
2243
2310
}
2244
2311
 
2245
2312
sub _write_sf {
2246
2313
    my $self = shift;
2247
2314
    my $sf = shift;
2248
 
    printf "TYPE:%s\n", $sf->primary_tag;
 
2315
    printf STDERR "TYPE:%s\n", $sf->primary_tag;
2249
2316
    return;
2250
2317
}
2251
2318
 
2252
2319
sub _write_sf_detail {
2253
2320
    my $self = shift;
2254
2321
    my $sf = shift;
2255
 
    printf "TYPE:%s\n", $sf->primary_tag;
 
2322
    printf STDERR "TYPE:%s\n", $sf->primary_tag;
2256
2323
    my @locs = $sf->location->each_Location;
2257
 
    printf "  %s,%s [%s]\n", $_->start, $_->end, $_->strand foreach @locs;
 
2324
    printf STDERR "  %s,%s [%s]\n", $_->start, $_->end, $_->strand foreach @locs;
2258
2325
    return;
2259
2326
}
2260
2327
 
2262
2329
    my $self = shift;
2263
2330
    my @sfs = @{shift || []};
2264
2331
    my $indent = shift || 0;
2265
 
    foreach my $sf (@sfs) {
2266
 
        my $label = '?';
2267
 
        if ($sf->has_tag('product')) {
2268
 
            ($label) = $sf->get_tag_values('product');
2269
 
        }
2270
 
        printf "%s%s $label\n", '  ' x $indent, $sf->primary_tag;
2271
 
        my @sub_sfs = $sf->sub_SeqFeature;
2272
 
        $self->_write_hier(\@sub_sfs, $indent+1);
 
2332
    if( $self->verbose > 0 ) {
 
2333
        foreach my $sf (@sfs) {
 
2334
            my $label = '?';
 
2335
            if ($sf->has_tag('product')) {
 
2336
                ($label) = $sf->get_tag_values('product');
 
2337
            }
 
2338
            printf STDERR "%s%s $label\n", '  ' x $indent, $sf->primary_tag;
 
2339
            my @sub_sfs = $sf->sub_SeqFeature;
 
2340
            $self->_write_hier(\@sub_sfs, $indent+1);
 
2341
        }
2273
2342
    }
2274
2343
}
2275
2344
 
2301
2370
               $inside = 0;
2302
2371
           }
2303
2372
       }
2304
 
       if ($self->verbose) {
2305
 
           print "    Checking containment:[$inside] (@container_coords) IN ($splice_uniq_str)\n";
 
2373
       if ($self->verbose > 0) {
 
2374
           printf STDERR "    Checking containment:[$inside] (@container_coords) IN ($splice_uniq_str)\n";
2306
2375
       }
2307
2376
       if ($inside) {
2308
2377
           # SCORE: matching (ss-scoords+2)/(n-container-ss-coords+2)
2390
2459
       $sf->isa("Bio::SeqFeatureI") || $self->throw("$sf NOT A SeqFeatureI");
2391
2460
       $sf->isa("Bio::FeatureHolderI") || $self->throw("$sf NOT A FeatureHolderI");
2392
2461
 
2393
 
       # so far, we only infer exons from mRNA
2394
2462
       my $type = $sf->primary_tag;
2395
2463
       next unless $type eq 'mRNA' or $type =~ /RNA/;
2396
2464
 
2397
2465
       # an mRNA from genbank will have a discontinuous location,
2398
2466
       # with each sub-location being equivalent to an exon
2399
2467
       my @locs = $sf->location;
 
2468
 
2400
2469
       if ($sf->location->isa("Bio::Location::SplitLocationI")) {
2401
2470
           @locs = $sf->location->each_Location;
2402
2471
       }
2403
 
       
 
2472
 
 
2473
       if (!@locs) {
 
2474
           use Data::Dumper;
 
2475
           print Dumper $sf;
 
2476
           $self->throw("ASSERTION ERROR: sf has no location objects");
 
2477
       }
 
2478
 
2404
2479
       # make exons from locations
2405
2480
       my @subsfs =
2406
2481
         map {
2407
2482
             my $subsf = Bio::SeqFeature::Generic->new(-location=>$_,
2408
2483
                                                       -primary_tag=>'exon');
 
2484
             ## Provide seq_id to new feature:
 
2485
             $subsf->seq_id($sf->seq_id) if $sf->seq_id;
 
2486
             $subsf->source_tag($sf->source_tag) if $sf->source_tag;
 
2487
             ## Transfer /locus_tag and /gene tag values to inferred
 
2488
             ## features.  TODO: Perhaps? this should not be done
 
2489
             ## indiscriminantly but rather by virtue of the setting
 
2490
             ## of group_tag.
 
2491
             foreach my $tag (grep /gene|locus_tag/, $sf->get_all_tags) {
 
2492
                 my @vals = $sf->get_tag_values($tag);
 
2493
                 $subsf->add_tag_value($tag, @vals);
 
2494
             }
 
2495
 
2409
2496
             my $locstr = 'exon::'.$self->_locstr($subsf);
2410
2497
 
2411
2498
             # re-use feature if type and location the same
2419
2506
         } @locs;
2420
2507
       
2421
2508
       # PARANOID CHECK
2422
 
       my $ok =
2423
 
         $self->_check_order_is_consistent(@subsfs);
2424
 
       if (!$ok) {
2425
 
           print "Unordered features:\n";
2426
 
           $self->_write_sf_detail($_) foreach @subsfs;
2427
 
           $self->throw("ASSERTION ERROR: inconsistent order");
2428
 
       }
 
2509
       $self->_check_order_is_consistent($sf->location->strand,@subsfs);
2429
2510
       #----
2430
2511
 
2431
2512
       $sf->location(Bio::Location::Simple->new());
 
2513
 
 
2514
       # we allow the exons to define the boundaries of the transcript
2432
2515
       $sf->add_SeqFeature($_, 'EXPAND') foreach @subsfs;
2433
2516
 
2434
2517
 
2514
2597
 
2515
2598
       $sf->isa("Bio::SeqFeatureI") || $self->throw("$sf NOT A SeqFeatureI");
2516
2599
       $sf->isa("Bio::FeatureHolderI") || $self->throw("$sf NOT A FeatureHolderI");
 
2600
       if ($self->verbose > 0) {
 
2601
           printf STDERR "    Checking $sf %s\n", $sf->primary_tag;
 
2602
       }
2517
2603
       
2518
2604
       if ($sf->primary_tag eq 'mRNA') {
2519
2605
           $self->problem(2,
2526
2612
           my @mrnas = ();
2527
2613
           foreach my $cds (@cdsl) {
2528
2614
               
2529
 
               my $ok;
2530
 
               $ok =
2531
 
                 $self->_check_order_is_consistent($cds->location->each_Location);
2532
 
               if (!$ok) {
2533
 
                   $self->_write_sf_detail($cds);
2534
 
                   $self->throw("inconsistent order");
2535
 
               }
2536
 
 
 
2615
               if ($self->verbose > 0) {
 
2616
                   print "    Inferring mRNA from CDS $cds\n";
 
2617
               }
 
2618
               $self->_check_order_is_consistent($cds->location->strand,$cds->location->each_Location);
 
2619
               
2537
2620
               my $loc = Bio::Location::Split->new;
2538
2621
               foreach my $cdsexonloc ($cds->location->each_Location) {
2539
2622
                   my $subloc =
2547
2630
                 Bio::SeqFeature::Generic->new(-location=>$loc,
2548
2631
                                               -primary_tag=>'mRNA');
2549
2632
               
2550
 
               $ok =
2551
 
                 $self->_check_order_is_consistent($mrna->location->each_Location);
2552
 
               if (!$ok) {
2553
 
                   $self->throw("inconsistent order");
2554
 
               }
 
2633
               ## Provide seq_id to new feature:
 
2634
               $mrna->seq_id($cds->seq_id) if $cds->seq_id;
 
2635
               $mrna->source_tag($cds->source_tag) if $cds->source_tag;
 
2636
 
 
2637
               $self->_check_order_is_consistent($mrna->location->strand,$mrna->location->each_Location);
 
2638
 
 
2639
               # make the mRNA hold the CDS; no EXPAND option,
 
2640
               # the CDS cannot be wider than the mRNA
2555
2641
               $mrna->add_SeqFeature($cds);
 
2642
 
2556
2643
               # mRNA steals children of CDS
2557
2644
               foreach my $subsf ($cds->get_SeqFeatures) {
2558
2645
                   $mrna->add_SeqFeature($subsf);
2609
2696
}
2610
2697
 
2611
2698
 
 
2699
# _check_order_is_consistent($strand,$ranges) RETURNS BOOL
 
2700
#
 
2701
# note: the value of this test is moot - there are many valid,
 
2702
# if unusual cases where it would flag an anomaly. for example
 
2703
# transpliced genes such as mod(mdg4) in dmel on AE003744, and
 
2704
# the following spliced gene on NC_001284:
 
2705
#
 
2706
#     mRNA            complement(join(20571..20717,21692..22086,190740..190761,
 
2707
#                     140724..141939,142769..142998))
 
2708
#                     /gene="nad5"
 
2709
#                     /note="trans-splicing, RNA editing"
 
2710
#                     /db_xref="GeneID:814567"
 
2711
#
 
2712
# note how the exons are not in order
 
2713
#  this will flag a level-3 warning, the user of this module
 
2714
#  can ignore this and deal appropriately with the resulting
 
2715
#  unordered exons
2612
2716
sub _check_order_is_consistent {
2613
2717
    my $self = shift;
 
2718
 
 
2719
    my $parent_strand = shift; # this does nothing..?
2614
2720
    my @ranges = @_;
2615
2721
    return unless @ranges;
 
2722
    my $rangestr =
 
2723
      join(" ",map{sprintf("[%s,%s]",$_->start,$_->end)} @ranges);
2616
2724
    my $strand = $ranges[0]->strand;
2617
2725
    for (my $i=1; $i<@ranges;$i++) {
2618
2726
        if ($ranges[$i]->strand != $strand) {
2619
 
            return 1; # mixed ranges - autopass
 
2727
            $self->problem(1,"inconsistent strands. Trans-spliced gene? Range: $rangestr");
 
2728
            return 1; 
 
2729
            # mixed ranges - autopass
 
2730
            # some mRNAs have exons on both strands; for
 
2731
            # example, the dmel mod(mdg4) gene which is
 
2732
            # trans-spliced (in actual fact two mRNAs)
2620
2733
        }
2621
2734
    }
 
2735
    my $pass = 1;
2622
2736
    for (my $i=1; $i<@ranges;$i++) {
2623
2737
        my $rangeP = $ranges[$i-1];
2624
2738
        my $range = $ranges[$i];
2625
 
#       if ($strand < 0) {
2626
 
#           if ($rangeP->end < $range->start) {
2627
 
#               return 0;
2628
 
#           }
2629
 
#       }
2630
 
#       else {
2631
2739
            if ($rangeP->start > $range->end) {
2632
 
                return 0;
2633
 
            }
2634
 
#       }
 
2740
                # failed - but still get one more chance..
 
2741
                $pass = 0;
 
2742
                $self->problem(2,"Ranges not in correct order. Strange ensembl genbank entry? Range: $rangestr");
 
2743
                last;
 
2744
            }
 
2745
    }
 
2746
    
 
2747
    if (!$pass) {
 
2748
        # sometimes (eg ensembl flavour genbank files)
 
2749
        # exons on reverse strand listed in reverse order
 
2750
        # eg join(complement(R1),...,complement(Rn))
 
2751
        # where R1 > R2
 
2752
        for (my $i=1; $i<@ranges;$i++) {
 
2753
            my $rangeP = $ranges[$i-1];
 
2754
            my $range = $ranges[$i];
 
2755
            if ($rangeP->end < $range->start) {
 
2756
                $self->problem(3,"inconsistent order. Range: $rangestr");
 
2757
                return 0;
 
2758
            }
 
2759
        }
2635
2760
    }
2636
2761
    return 1; # pass
2637
2762
}
2670
2795
    
2671
2796
}
2672
2797
 
 
2798
1;