Partitioning 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
                     /gene="noc"
                     /locus_tag="CG4491"
                     /note="last curated on Thu Dec 13 16:51:32 PST 2001"
                     /map="35B2-35B2"
                     /db_xref="FLYBASE:FBgn0005771"
     mRNA            join(20111..20584,20887..23268)
                     /gene="noc"
                     /locus_tag="CG4491"
                     /product="CG4491-RA"
                     /db_xref="FLYBASE:FBgn0005771"

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 disconnected.

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

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:

  gene
    mRNA
      CDS
    mRNA
      CDS
    mRNA
      CDS

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

  gene
    mRNA
      CDS

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

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

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

  $unflattener->feature_from_splitloc(-seq=>$seq);

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

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

  gene
    CDS

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

  gene
    mRNA
      CDS

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

   $unflattener->infer_mRNA_from_CDS(-seq=>$seq);

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

Customising the grouping of features

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

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)
                     /locus_tag="CG32954"
                     /note="CG32954 gene product from transcript CG32954-RA"
                                                 ^^^^^^^^^^^^^^^^^^^^^^^^^^^
                     /codon_start=1
                     /product="CG32954-PA"
                     /protein_id="AAF53403.1"
                     /db_xref="GI:7298167"
                     /db_xref="FLYBASE:FBgn0052954"
                     /translation="MSFTLTNKNVIFVAGLGGIGLDTSKELLKRDLKNLVILDRIENP..."

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

  /product="CG32954-RA"

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 this.

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.

       $unflattener->unflatten_seq(-seq=>$seq,
                                 -group_tag=>'locus_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) {
                                             $self->throw("UNRESOLVABLE");
                                         }
                                         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) {
                                             $self->throw("UNRESOLVABLE");
                                         }
                                         elsif (@container_sfs == 1) {
                                             # we got it!
                                             return $container_sfs[0];
                                         }
                                         else {
                                             $self->throw("AMBIGUOUS");
                                         }
                                     }
                                 });

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

Parsing mRNA records

NOT YET DONE - IN PROGRESS!!!

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...?

Mailing Lists

  bioperl-l@bioperl.org                         - General discussion
  http://bioperl.org/wiki/Mailing_lists  - About the mailing lists

Support

bioperl-l@bioperl.org

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

  https://github.com/bioperl/bioperl-live/issues

new

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

Arguments

  -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)

seq

 Title   : seq
 Usage   : $unflattener->seq($newval)
 Function: 
 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()

group_tag

 Title   : group_tag
 Usage   : $unflattener->group_tag($newval)
 Function: 
 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
                     /gene="noc"
                     /locus_tag="CG4491"
                     /note="last curated on Thu Dec 13 16:51:32 PST 2001"
                     /map="35B2-35B2"
                     /db_xref="FLYBASE:FBgn0005771"
     mRNA            join(20111..20584,20887..23268)
                     /gene="noc"
                     /locus_tag="CG4491"
                     /product="CG4491-RA"
                     /db_xref="FLYBASE:FBgn0005771"

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

partonomy

 Title   : partonomy
 Usage   : $unflattener->partonomy({mRNA=>'gene', CDS=>'mRNA')
 Function: 
 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.

structure_type

 Title   : structure_type
 Usage   : $unflattener->structure_type($newval)
 Function: 
 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 record.

Type 0 (DEFAULT)

typically contains

  source
  gene
  mRNA
  CDS
    

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

  gene
    mRNA
    CDS
      exon
    

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

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

  $unflattener->feature_from_splitloc(-seq=>$seq);
    

Type 1

typically contains

  source
  gene
  CDS
  exon [optional]
  intron [optional]
    

there are no mRNA features

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

  gene
    CDS
      exon
      intron
    

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

get_problems

 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

clear_problems

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

report_problems

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

ignore_problems

 Title   : ignore_problems
 Usage   : $obj->ignore_problems();
 Function:
 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.

error_threshold

 Title   : error_threshold
 Usage   : $obj->error_threshold($severity)
 Function: 
 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.

unflatten_seq

 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 explanation.

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

Arguments

  -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'

unflatten_groups

 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

Arguments

  -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.

unflatten_group

 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

Arguments

  -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.

feature_from_splitloc

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

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

Arguments:

  -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

infer_mRNA_from_CDS

 Title   : infer_mRNA_from_CDS
 Usage   :
 Function:
 Example :
 Returns : 
 Args    :

given a "type 1" containment hierarchy

  gene
    CDS
      exon

this will infer the uniform "type 0" containment hierarchy

  gene
    mRNA
      CDS
      exon

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 ftp.ncbi.nih.gov/genomes/Schizosaccharomyces_pombe/)

remove_types

 Title   : remove_types
 Usage   : $unf->remove_types(-seq=>$seq, -types=>["mRNA"]);
 Function:
 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 ftp.ncbi.nih.gov/genomes/Schizosaccharomyces_pombe/ UNLESS the bogus mRNAs in these records are removed (or changed to a different type) - they just confuse things too much