Provided by: libbio-perl-perl_1.7.2-2_all bug


       Bio::SeqFeature::Tools::Unflattener - turns flat list of genbank-sourced features into a
       nested SeqFeatureI hierarchy


         # standard / generic use - unflatten a genbank record
         use Bio::SeqIO;
         use Bio::SeqFeature::Tools::Unflattener;

         # generate an Unflattener object
         $unflattener = Bio::SeqFeature::Tools::Unflattener->new;

         # first fetch a genbank SeqI object
         $seqio =
         my $out =
         while ($seq = $seqio->next_seq()) {

           # get top level unflattended SeqFeatureI objects

           @top_sfs = $seq->get_SeqFeatures;
           foreach my $sf (@top_sfs) {
               # do something with top-level features (eg genes)


       Most GenBank entries for annotated genomic DNA contain a flat list of features. These
       features can be parsed into an equivalent flat list of Bio::SeqFeatureI objects using the
       standard Bio::SeqIO classes. However, it is often desirable to unflatten this list into
       something resembling actual gene models, in which genes, mRNAs and CDSs are nested
       according to the nature of the gene model.

       The BioPerl object model allows us to store these kind of associations between SeqFeatures
       in containment hierarchies -- any SeqFeatureI object can contain nested SeqFeatureI
       objects. The Bio::SeqFeature::Tools::Unflattener object facilitates construction of these
       hierarchies from the underlying GenBank flat-feature-list representation.

       For example, if you were to look at a typical GenBank DNA entry, say, AE003644, you would
       see a flat list of features:


         gene CG4491
         mRNA CG4491-RA
         CDS CG4491-PA

         gene tRNA-Pro
         tRNA tRNA-Pro

         gene CG32954
         mRNA CG32954-RA
         mRNA CG32954-RC
         mRNA CG32954-RB
         CDS CG32954-PA
         CDS CG32954-PB
         CDS CG32954-PC

       These features have sequence locations, but it is not immediately clear how to write code
       such that each mRNA is linked to the appropriate CDS (other than relying on IDs which is
       very bad)

       We would like to convert the above list into the containment hierarchy, shown below:

           mRNA CG4491-RA
             CDS CG4491-PA
           tRNA tRNA-Pro
           mRNA CG32954-RA
             CDS CG32954-PA
           mRNA CG32954-RC
             CDS CG32954-PC
           mRNA CG32954-RB
             CDS CG32954-PB

       Where each feature is nested underneath its container. Note that exons have been
       automatically inferred (even for tRNA genes).

       We do this using a call on a Bio::SeqFeature::Tools::Unflattener object

         @sfs = $unflattener->unflatten_seq(-seq=>$seq);

       This would return a list of the top level (i.e. container) SeqFeatureI objects - in this
       case, genes. Other top level features are possible; for instance, the source feature which
       is always present, and other features such as variation or misc_feature types.

       The containment hierarchy can be accessed using the get_SeqFeature() call on any feature
       object - see Bio::SeqFeature::FeatureHolderI.  The following code will traverse the
       containment hierarchy for a feature:

         sub traverse {
           $sf = shift;   #  $sf isa Bio::SeqfeatureI

           # something with $sf!

           # depth first traversal of containment tree
           @contained_sfs = $sf->get_SeqFeatures;
           traverse($_) foreach @contained_sfs;

       Once you have built the hierarchy, you can do neat stuff like turn the features into
       'rich' feature objects (eg Bio::SeqFeature::Gene::GeneStructure) or convert to a suitable
       format such as GFF3 or chadoxml (after mapping to the Sequence Ontology); this step is not
       described here.


       Due to the quixotic nature of how features are stored in GenBank/EMBL/DDBJ, there is no
       guarantee that the default behaviour of this module will produce perfect results.
       Sometimes it is hard or impossible to build a correct containment hierarchy if the
       information provided is simply too lossy, as is often the case. If you care deeply about
       your data, you should always manually inspect the resulting containment hierarchy; you may
       have to customise the algorithm for building the hierarchy, or even manually tweak the
       resulting hierarchy. This is explained in more detail further on in the document.

       However, if you are satisfied with the default behaviour, then you do not need to read any
       further. Just make sure you set the parameter use_magic - this will invoke incantations
       which will magically produce good results no matter what the idiosyncracies of the
       particular GenBank record in question.

       For example


       The success of this depends on the phase of the moon at the time the entry was submitted
       to GenBank. Note that the magical recipe is being constantly improved, so the results of
       invoking magic may vary depending on the bioperl release.

       If you are skeptical of magic, or you wish to exact fine grained control over how the
       entry is unflattened, or you simply wish to understand more about how this crazy stuff
       works, then read on!


       Occasionally the Unflattener will have problems with certain records. For example, the
       record may contain inconsistent data - maybe there is an exon entry that has no
       corresponding mRNA location.

       The default behaviour is to throw an exception reporting the problem, if the problem is
       relatively serious - for example, inconsistent data.

       You can exert more fine grained control over this - perhaps you want the Unflattener to do
       the best it can, and report any problems. This can be done - refer to the methods.






       This is the default algorithm; you should be able to override any part of it to customise.

       The core of the algorithm is in two parts

       Partitioning the flat feature list into groups
       Resolving the feature containment hierarchy for each group

       There are other optional steps after the completion of these two steps, such as inferring
       exons; we now describe in more detail what is going on.

   Partitioning into groups
       First of all the flat feature list is partitioned into groups.

       The default way of doing this is to use the gene attribute; if we look at two features
       from GenBank accession AE003644.3:

            gene            20111..23268
                            /note="last curated on Thu Dec 13 16:51:32 PST 2001"
            mRNA            join(20111..20584,20887..23268)

       Both these features share the same /gene tag which is "noc", so they correspond to the
       same gene model (the CDS feature is not shown, but this also has a tag-value /gene="noc").

       Not all groups need to correspond to gene models, but this is the most common use case;
       later on we shall describe how to customise the grouping.

       Sometimes other tags have to be used; for instance, if you look at the entire record for
       AE003644.3 you will see you actually need the use the /locus_tag attribute. This attribute
       is actually not present in most records!

       You can override this:

         $collection->unflatten_seq(-seq=>$seq, -group_tag=>'locus_tag');

       Alternatively, if you -use_magic, the object will try and make a guess as to what the
       correct group_tag should be.

       At the end of this step, we should have a list of groups - there is no structure within a
       group; the group just serves to partition the flat features. For the example data above,
       we would have the following groups.

         [ source ]
         [ gene mRNA CDS ]
         [ gene mRNA CDS ]
         [ gene mRNA CDS ]
         [ gene mRNA mRNA mRNA CDS CDS CDS ]

       Multicopy Genes

       Multicopy genes are usually rRNAs or tRNAs that are duplicated across the genome. Because
       they are functionally equivalent, and usually have the same sequence, they usually have
       the same group_tag (ie gene symbol); they often have a /note tag giving copy number. This
       means they will end up in the same group. This is undesirable, because they are spatially

       There is another step, which involves splitting spatially disconnected groups into
       distinct groups

       this would turn this

        [gene-rrn3 rRNA-rrn3 gene-rrn3 rRNA-rrn3]

       into this

        [gene-rrn3 rRNA-rrn3] [gene-rrn3 rRNA-rrn3]

       based on the coordinates

       What next?

       The next step is to add some structure to each group, by making containment hierarchies,
       trees that represent how the features interrelate

   Resolving the containment hierarchy
       After the grouping is done, we end up with a list of groups which probably contain
       features of type 'gene', 'mRNA', 'CDS' and so on.

       Singleton groups (eg the 'source' feature) are ignored at this stage.

       Each group is itself flat; we need to add an extra level of organisation. Usually this is
       because different spliceforms (represented by the 'mRNA' feature) can give rise to
       different protein products (indicated by the 'CDS' feature). We want to correctly
       associate mRNAs to CDSs.

       We want to go from a group like this:

         [ gene mRNA mRNA mRNA CDS CDS CDS ]

       to a containment hierarchy like this:


       In which each CDS is nested underneath the correct corresponding mRNA.

       For entries that contain no alternate splicing, this is simple; we know that the group

         [ gene mRNA CDS ]

       Must resolve to the tree


       How can we do this in entries with alternate splicing? The bad news is that there is no
       guaranteed way of doing this correctly for any GenBank entry. Occasionally the submission
       will have been done in such a way as to reconstruct the containment hierarchy. However,
       this is not consistent across databank entries, so no generic solution can be provided by
       this object. This module does provide the framework within which you can customise a
       solution for the particular dataset you are interested in - see later.

       The good news is that there is an inference we can do that should produce pretty good
       results the vast majority of the time. It uses splice coordinate data - this is the
       default behaviour of this module, and is described in detail below.

   Using splice site coordinates to infer containment
       If an mRNA is to be the container for a CDS, then the splice site coordinates (or intron
       coordinates, depending on how you look at it) of the CDS must fit inside the splice site
       coordinates of the mRNA.

       Ambiguities can still arise, but the results produced should still be reasonable and
       consistent at the sequence level. Look at this fake example:

         mRNA    XXX---XX--XXXXXX--XXXX         join(1..3,7..8,11..16,19..23)
         mRNA    XXX-------XXXXXX--XXXX         join(1..3,11..16,19..23)
         CDS                 XXXX--XX           join(13..16,19..20)
         CDS                 XXXX--XX           join(13..16,19..20)

       [obviously the positions have been scaled down]

       We cannot unambiguously match mRNA with CDS based on splice sites, since both CDS share
       the splice site locations 16^17 and 18^19. However, the consequences of making a wrong
       match are probably not very severe. Any annotation data attached to the first CDS is
       probably identical to the seconds CDS, other than identifiers.

       The default behaviour of this module is to make an arbitrary call where it is ambiguous
       (the mapping will always be bijective; i.e. one mRNA -> one CDS).

       [TODO: NOTE: not tested on EMBL data, which may not be bijective; ie two mRNAs can share
       the same CDS??]

       This completes the building of the containment hierarchy; other optional step follow


   Inferring exons from mRNAs
       This step always occurs if -use_magic is invoked.

       In a typical GenBank entry, the exons are implicit. That is they can be inferred from the
       mRNA location.

       For example:

            mRNA            join(20111..20584,20887..23268)

       This tells us that this particular transcript has two exons. In bioperl, the mRNA feature
       will have a 'split location'.

       If we call


       This will generate the necessary exon features, and nest them under the appropriate mRNAs.
       Note that the mRNAs will no longer have split locations - they will have simple locations
       spanning the extent of the exons. This is intentional, to avoid redundancy.

       Occasionally a GenBank entry will have both implicit exons (from the mRNA location) and
       explicit exon features.

       In this case, exons will still be transferred. Tag-value data from the explicit exon will
       be transferred to the implicit exon. If exons are shared between mRNAs these will be
       represented by different objects. Any inconsistencies between implicit and explicit will
       be reported.

       tRNAs and other noncoding RNAs

       exons will also be generated from these features

   Inferring mRNAs from CDS
       Some GenBank entries represent gene models using features of type gene, mRNA and CDS; some
       entries just use gene and CDS.

       If we only have gene and CDS, then the containment hierarchies will look like this:


       If we want the containment hierarchies to be uniform, like this


       Then we must create an mRNA feature. This will have identical coordinates to the CDS. The
       assumption is that there is either no untranslated region, or it is unknown.

       To do this, we can call


       This is taken care of automatically, if -use_magic is invoked.


   Customising the grouping of features
       The default behaviour is suited mostly to building models of protein coding genes and
       noncoding genes from genbank genomic DNA submissions.

       You can change the tag used to partition the feature by passing in a different group_tag
       argument - see the unflatten_seq() method

       Other behaviour may be desirable. For example, even though SNPs (features of type
       'variation' in GenBank) are not actually part of the gene model, it may be desirable to
       group SNPs that overlap or are nearby gene models.

       It should certainly be possible to extend this module to do this. However, I have yet to
       code this part!!! If anyone would find this useful let me know.

       In the meantime, you could write your own grouping subroutine, and feed the results into
       unflatten_groups() [see the method documentation below]

   Customising the resolution of the containment hierarchy
       Once the flat list of features has been partitioned into groups, the method
       unflatten_group() is called on each group to build a tree.

       The algorithm for doing this is described above; ambiguities are resolved by using splice
       coordinates. As discussed, this can be ambiguous.

       Some submissions may contain information in tags/attributes that hint as to the mapping
       that needs to be made between the features.

       For example, with the Drosophila Melanogaster release 3 submission, we see that CDS
       features in alternately spliced mRNAs have a form like this:

            CDS             join(145588..145686,145752..146156,146227..146493)
                            /note="CG32954 gene product from transcript CG32954-RA"

       Here the /note tag provides the clue we need to link CDS to mRNA (highlighted with ^^^^).
       We just need to find the mRNA with the tag


       I have no idea how consistent this practice is across submissions; it is consistent for
       the fruitfly genome submission.

       We can customise the behaviour of unflatten_group() by providing our own resolver method.
       This obviously requires a bit of extra programming, but there is no way to get around

       Here is an example of how to pass in your own resolver; this example basically checks the
       parent (container) /product tag to see if it matches the required string in the child
       (contained) /note tag.

                                        -resolver_method=>sub {
                                            my $self = shift;
                                            my ($sf, @candidate_container_sfs) = @_;
                                            if ($sf->has_tag('note')) {
                                                my @notes = $sf->get_tag_values('note');
                                                my @trnames = map {/from transcript\s+(.*)/;
                                                                   $1} @notes;
                                                @trnames = grep {$_} @trnames;
                                                my $trname;
                                                if (@trnames == 0) {
                                                elsif (@trnames == 1) {
                                                    $trname = $trnames[0];
                                                else {
                                                    $self->throw("AMBIGUOUS: @trnames");
                                                my @container_sfs =
                                                  grep {
                                                      my ($product) =
                                                        $_->has_tag('product') ?
                                                          $_->get_tag_values('product') :
                                                      $product eq $trname;
                                                  } @candidate_container_sfs;
                                                if (@container_sfs == 0) {
                                                elsif (@container_sfs == 1) {
                                                    # we got it!
                                                    return $container_sfs[0];
                                                else {

       the resolver method is only called when there is more than one spliceform.

   Parsing mRNA records
       Some of the entries in sequence databanks are for mRNA sequences as well as genomic DNA.
       We may want to build models from these too.


       Open question - what would these look like?

       Ideally we would like a way of combining a mRNA record with the corresponding SeFeature
       entry from the appropriate genomic DNA record. This could be problemmatic in some cases -
       for example, the mRNA sequences may not match 100% (due to differences in strain, assembly
       problems, sequencing problems, etc). What then...?


       Feature table description


   Mailing Lists
       User feedback is an integral part of the evolution of this and other Bioperl modules. Send
       your comments and suggestions preferably to the Bioperl mailing lists  Your participation
       is much appreciated.
                         - General discussion  - About the mailing lists

       Please direct usage questions or support issues to the mailing list:

       rather than to the module maintainer directly. Many experienced and reponsive experts will
       be able look at the problem and quickly address it. Please include a thorough description
       of the problem with code and data examples if at all possible.

   Reporting Bugs
       report bugs to the Bioperl bug tracking system to help us keep track the bugs and their
       resolution.  Bug reports can be submitted via the web:

AUTHOR - Chris Mungall



       The rest of the documentation details each of the object methods. Internal methods are
       usually preceded with a _

        Title   : new
        Usage   : $unflattener = Bio::SeqFeature::Tools::Unflattener->new();
        Function: constructor
        Example :
        Returns : a new Bio::SeqFeature::Tools::Unflattener
        Args    : see below


         -seq       : A L<Bio::SeqI> object (optional)
                      the sequence to unflatten; this can also be passed in
                      when we call unflatten_seq()

         -group_tag : a string representing the /tag used to partition flat features
                      (see discussion above)

        Title   : seq
        Usage   : $unflattener->seq($newval)
        Example :
        Returns : value of seq (a Bio::SeqI)
        Args    : on set, new value (a Bio::SeqI, optional)

       The Bio::SeqI object should hold a flat list of Bio::SeqFeatureI objects; this is the list
       that will be unflattened.

       The sequence object can also be set when we call unflatten_seq()

        Title   : group_tag
        Usage   : $unflattener->group_tag($newval)
        Example :
        Returns : value of group_tag (a scalar)
        Args    : on set, new value (a scalar or undef, optional)

       This is the tag that will be used to collect elements from the flat feature list into
       groups; for instance, if we look at two typical GenBank features:

            gene            20111..23268
                            /note="last curated on Thu Dec 13 16:51:32 PST 2001"
            mRNA            join(20111..20584,20887..23268)

       We can see that these comprise the same gene model because they share the same /gene
       attribute; we want to collect these together in groups.

       Setting group_tag is optional. The default is to use 'gene'. In the example above, we
       could also use /locus_tag

        Title   : partonomy
        Usage   : $unflattener->partonomy({mRNA=>'gene', CDS=>'mRNA')
        Example :
        Returns : value of partonomy (a scalar)
        Args    : on set, new value (a scalar or undef, optional)

       A hash representing the containment structure that the seq_feature nesting should conform
       to; each key represents the contained (child) type; each value represents the container
       (parent) type.

        Title   : structure_type
        Usage   : $unflattener->structure_type($newval)
        Example :
        Returns : value of structure_type (a scalar)
        Args    : on set, new value (an int or undef, optional)

       GenBank entries conform to different flavours, or structure types. Some have mRNAs, some
       do not.

       Right now there are only two base structure types defined. If you set the structure type,
       then appropriate unflattening action will be taken.  The presence or absence of explicit
       exons does not affect the structure type.

       If you invoke -use_magic then this will be set automatically, based on the content of the

       Type 0 (DEFAULT)
           typically contains


           with this structure type, we want the seq_features to be nested like this


           exons and introns are implicit from the mRNA 'join' location

           to get exons from the mRNAs, you will need this call (see below)


       Type 1
           typically contains

             exon [optional]
             intron [optional]

           there are no mRNA features

           with this structure type, we want the seq_features to be nested like this


           exon and intron may or may not be present; they may be implicit from the CDS 'join'

        Title   : get_problems
        Usage   : @probs = get_problems()
        Function: Get the list of problem(s) for this object.
        Example :
        Returns : An array of [severity, description] pairs
        Args    :

       In the course of unflattening a record, problems may occur. Some of these problems are
       non-fatal, and can be ignored.

       Problems are represented as arrayrefs containing a pair [severity, description]

       severity is a number, the higher, the more severe the problem

       the description is a text string

        Title   : clear_problems
        Usage   :
        Function: resets the problem list to empty
        Example :
        Returns :
        Args    :

        Title   : report_problems
        Usage   : $unflattener->report_problems(\*STDERR);
        Example :
        Returns :
        Args    : FileHandle (defaults to STDERR)

        Title   : ignore_problems
        Usage   : $obj->ignore_problems();
        Example :
        Returns :
        Args    :

       Unflattener is very particular about problems it finds along the way. If you have set the
       error_threshold such that less severe problems do not cause exceptions, Unflattener still
       expects you to report_problems() at the end, so that the user of the module is aware of
       any inconsistencies or problems with the data. In fact, a warning will be produced if
       there are unreported problems. To silence, this warning, call the ignore_problems() method
       before the Unflattener object is destroyed.

        Title   : error_threshold
        Usage   : $obj->error_threshold($severity)
        Example :
        Returns : value of error_threshold (a scalar)
        Args    : on set, new value (an integer)

       Sets the threshold above which errors cause this module to throw an exception. The default
       is 0; all problems with a severity > 0 will cause an exception.

       If you raise the threshold to 1, then the unflattening process will be more lax; problems
       of severity==1 are generally non-fatal, but may indicate that the results should be
       inspected, for example, to make sure there is no data loss.

        Title   : unflatten_seq
        Usage   : @sfs = $unflattener->unflatten_seq($seq);
        Function: turns a flat list of features into a list of holder features
        Example :
        Returns : list of Bio::SeqFeatureI objects
        Args    : see below

       partitions a list of features then arranges them in a nested tree; see above for full

       note - the Bio::SeqI object passed in will be modified


         -seq   :          a Bio::SeqI object; must contain Bio::SeqFeatureI objects
                           (this is optional if seq has already been set)

         -use_magic:       if TRUE (ie non-zero) then magic will be invoked;
                           see discussion above.

         -resolver_method: a CODE reference
                           see the documentation above for an example of
                           a subroutine that can be used to resolve hierarchies
                           within groups.

                           this is optional - if nothing is supplied, a default
                           subroutine will be used (see below)

         -group_tag:       a string
                           [ see the group_tag() method ]
                           this overrides the default group_tag which is 'gene'

        Title   : unflatten_groups
        Usage   :
        Function: iterates over groups, calling unflatten_group() [see below]
        Example :
        Returns : list of Bio::SeqFeatureI objects that are holders
        Args    : see below


         -groups:          list of list references; inner list is of Bio::SeqFeatureI objects
                           e.g.  ( [$sf1], [$sf2, $sf3, $sf4], [$sf5, ...], ...)

         -resolver_method: a CODE reference
                           see the documentation above for an example of
                           a subroutine that can be used to resolve hierarchies
                           within groups.

                           this is optional - a default subroutine will be used

       NOTE: You should not need to call this method, unless you want fine grained control over
       how the unflattening process.

        Title   : unflatten_group
        Usage   :
        Function: nests a group of features into a feature containment hierarchy
        Example :
        Returns : Bio::SeqFeatureI objects that holds other features
        Args    : see below


         -group:           reference to list of Bio::SeqFeatureI objects

         -resolver_method: a CODE reference
                           see the documentation above for an example of
                           a subroutine that can be used to resolve hierarchies
                           within groups

                           this is optional - a default subroutine will be used

       NOTE: You should not need to call this method, unless you want fine grained control over
       how the unflattening process.

        Title   : feature_from_splitloc
        Usage   : $unflattener->feature_from_splitloc(-features=>$sfs);
        Example :
        Returns :
        Args    : see below

       At this time all this method does is generate exons for mRNA or other RNA features


         -feature:    a Bio::SeqFeatureI object (that conforms to Bio::FeatureHolderI)
         -seq:        a Bio::SeqI object that contains Bio::SeqFeatureI objects
         -features:   an arrayref of Bio::SeqFeatureI object

        Title   : infer_mRNA_from_CDS
        Usage   :
        Example :
        Returns :
        Args    :

       given a "type 1" containment hierarchy


       this will infer the uniform "type 0" containment hierarchy


       all the children of the CDS will be moved to the mRNA

       a "type 2" containment hierarchy is mixed type "0" and "1" (for example, see

        Title   : remove_types
        Usage   : $unf->remove_types(-seq=>$seq, -types=>["mRNA"]);
        Example :
        Returns :
        Args    :

       removes features of a set type

       useful for pre-filtering a genbank record; eg to get rid of STSs

       also, there is no way to unflatten
       UNLESS the bogus mRNAs in these records are removed (or changed to a different type) -
       they just confuse things too much